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ABSTRACT 

The  objective  of  Bayesian  belief  propagation  in  this  paper  is  to  perform  an  interactive  status  assessment  of 
geosynchronous  satellites  as  each  new  data  point  for  the  photometric  brightness  becomes  available  during  the 
synoptic  search  performed  by  a  space-based  sensor  as  a  part  of  its  routine  metric  mission.  The  calculations  are 
performed  by  using  a  dimensionless  ratio  of  observed  photometric  brightness  to  its  predicted  brightness.  The 
brightness  predictions  can  be  obtained  using  any  analytical  model  chosen  by  the  user.  The  inference  for  a  level  of 
confidence  in  the  statistical  assessment  is  performed  on  the  basis  of  propagated  values  for  belief  within  a  cluster  of 
satellites  that  are  located  within  a  close  proximity  to  each  other.  This  is  meant  to  render  the  assessment  to  be  as 
independent  of  assumptions  and  algorithms  utilized  in  the  analytical  model  as  possible;  and  to  mitigate  the  effect  of 
bias  that  could  be  introduced  by  the  choice  of  analytical  model.  It  considers  that  a  model  performs  predictions  based 
on  the  geometry  of  observation  conditions  and  any  information  that  could  have  been  extracted  by  the  inversion  of 
prior  data  on  its  photometric  brightness.  Thus,  if  there  is  a  statistical  change  in  the  predictive  error  for  a  single 
satellite  or  a  pair  of  satellites,  while  remaining  unchanged  for  the  rest,  there  is  higher  likelihood  of  anomaly  in  either 
the  operational  status  of  that  satellite  or  an  error  in  object  correlation  (i.e.  cross-tag). 

The  algorithm  in  this  paper  uses  a  first  order  Markov  chain  model  to  compute  a  conditional  probability  value  for  the 
satellite  status  to  be  nominal  or  anomalous  (i.e.,  NOM  or  ANOM),  given  its  latest  photometry  observation.  This 
calculation  is  repeated  as  data  for  each  new  observation  becomes  available.  Also,  it  is  performed  for  each  satellite 
(member)  that  belongs  to  a  geosynchronous  cluster  (group).  This  provides  a  sequence  of  conditional  probability 
values  for  each  member  in  a  group.  This  information  is  mapped  to  a  tree-like  directed  graph  (i.e.  Bayesian  network) 
of  nodes.  There  is  a  node  for  each  new  data  point  and  it  represents  a  hypothesis  test  for  whether  the  member  status  is 
NOM  or  ANOM  at  the  time  of  that  observation.  The  conditional  probability  values  for  the  status  of  each  member  in 
the  group  are  utilized  in  order  to  compute  the  marginal  probability  (i.e.  belief)  that  the  status  of  an  individual 
member  is  NOM  or  ANOM.  The  propagation  of  belief  summarizes  all  preceding  computations  of  belief  in  order  to 
determine  a  level  of  confidence  in  the  statistical  assessment  for  its  use  by  an  analyst. 
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1.  INTRODUCTION 


The  method  reported  in  this  work  is  constructed  as  teamwork  of  three  analytical  procedures;  namely  (1)  A  procedure 
to  perform  inversion  of  non-resolved  photometry  data  for  the  brightness  of  three-axis  stabilized  geosynchronous 
satellites.  This  inversion  solves  for  the  individual  albedo-area  products  for  the  satellite  body  and  its  solar  panels  as  a 
function  of  their  respective  observation  geometry.  The  observation  geometry  is  defined  separately  for  the  body  and 
panel  in  terms  of  their  respective  angles  for  the  direction  of  illumination  and  the  sensor  line  of  sight.  This  procedure 
is  denoted  as  Inversion  Model;  (2)  A  procedure  to  predict  the  brightness  of  a  geosynchronous  satellite  at  a  given 
future  epoch  using  the  computed  values  for  the  albedo-area  products  for  the  satellite  body  and  its  solar  panels.  This 
procedure  is  denoted  as  Predictive  Model;  (3)  A  method  to  perform  statistical  assessment  in  near-real  time.  This  is 
denoted  as  Statistics  Model  and  it  is  the  focus  of  the  present  paper.  The  Inversion  and  Predictive  Models  are 
collectively  denoted  as  Photometry  Model.  The  notation  used  is  defined  in  Appendix  A  in  addition  to  the  first  time  it 
appears  in  the  technical  description. 

The  Statistics  Model  is  independent  of  the  Inversion  Model  and  the  Predictive  Model.  Each  model  is  constructed  as 
a  reusable  entity.  The  orchestration  of  teamwork  between  the  three  models  is  performed  within  the  Statistics  Model 
and  it  is  denoted  as  Control  Logic.  The  flow  chart  for  the  Control  Logic  is  shown  in  Section  16.  The  Photometry 
Model  is  assumed  to  be  imperfect  but  invariant.  It  is  imperfect  in  that  its  accuracy  changes  from  one  satellite  to 
another  and  from  one  observation  geometry  to  another.  It  is  invariant  in  that  its  analytical  competence  is  fixed.  This 
work  uses  the  Photometry  Model  reported  in  References  1-3.  It  is  constructed  using  the  principle  of  material  frame 
indifference,  which  allows  the  same  mathematics  to  be  useful  for  the  processing  of  ground-based  or  space-based 
sensor  data.  This  model  is  implemented  as  an  object-oriented  software  application  [Reference  25].  The  Statistics 
Model  is  newly  implemented  as  its  reusable  entity  within  the  same  software. 

The  present  work  is  constructed  such  that  its  models  are  useful  for  the  processing  of  brightness  data  collected  during 
the  routine  metrics  mission  performed  by  the  existing  ground  or  space-based  sensors  in  the  space  surveillance 
network.  This  is  a  synoptic  search  operation  and  it  is  the  workhorse  for  the  maintenance  of  geosynchronous  satellite 
catalog.  The  synoptic  search  data  includes  both  angles-only  metric  data  and  single-point  visual,  panchromatic 
brightness  data.  The  positional  information  for  the  satellites  is  derived  from  the  angles-only  data  and  is  utilized 
today  for  satellite  catalog  maintenance.  The  optical  cross-section  and  orientation  information  for  the  satellites  can  be 
derived  from  the  single  point  panchromatic  brightness  data,  but  this  information  is  neither  used  nor  calculated  at 
present  even  though  the  metric  data  is  not  sufficient  by  itself  to  maintain  accurate  satellite  catalog  due  to  the 
overcrowding  of  the  geosynchronous  (GEO)  belt  [e.g.  Reference  12].  There  are  two  reasons  why  the  single  point 
brightness  data  (or  just  ‘Brightness  data’  or  ‘Brightness’  for  short)  is  not  used  as  follows: 

The  first  reason  is  that  photometry  analyses  are  traditionally  performed  using  signature  data  (denoted  as  ‘Signature 
Data’  in  this  work)  collected  using  dedicated  ground-based  sensors.  Such  Signature  Data  is  typically  at  a  rate  of  one 
or  more  observations  per  minute.  It  extends  for  several  minutes  or  even  hours,  vividly  displaying  a  clear  temporal 
character  of  satellite  signatures  to  a  human  eye  [Reference  29].  On  the  contrary.  Brightness  Data  is  collected  at  a 
rate  two  orders  of  magnitude  slower,  hiding  the  same  knowledge  in  bits  and  pieces  in  the  sparseness.  Considering 
that  the  solar  phase  angle  changes  15°  per  hour,  if  each  observation  was  separated  by  90-minutes  (which  is  about 
equal  to  the  orbital  period  for  a  notional  space-based  sensor  in  LEO),  the  successive  points  in  the  Brightness  data 
may  be  considered  to  be  separated  by  about  22.5°  in  phase  angle.  In  addition,  there  is  no  data  during  daytime  (Ligure 
1.1).  A  new  formulation  of  mathematics  is  required  to  quilt  this  knowledge  together,  particularly  in  the  case  of 
space-based  data. 

The  second  reason  is  that  the  mathematics  of  inversion  of  photometry  comprises  more  unknowns  than  the  number  of 
independent  equations.  Thus,  a  mathematical  proof  cannot  be  offered  for  uniqueness  of  solution,  except  in  special 
situations.  This  lack  of  uniqueness  translates  into  lack  of  confidence  in  the  inference  derived  from  the  solution.  A 
new  formulation  of  mathematics  is  required  to  gain  the  requisite  level  of  confidence  through  a  mathematical 
evolution  of  belief  (i.e.  belief  propagation  or  BP)  derived  from  a  sequence  of  inversions  of  Brightness  Data. 

The  first  reason  can  be  addressed  by  formulating  the  governing  equations  for  inversion  of  Brightness  Data  using  the 
principle  of  material  frame  indifference  (PMLI)  [References  1-3].  This  work  draws  upon  the  hidden  richness  of 
information  in  the  Brightness  data  as  compared  to  the  traditional  Signature  Data.  Signature  Data  has  higher  frame 
rate  that  spans  the  full  range  of  longitudinal  phase  angles,  but  for  a  single  value  of  solar  declination.  Brightness  data 
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has  a  much  slower  frame  rate,  but  it  samples  the  full  range  of  longitudinal  phase  angle,  solar  declination  and,  in  the 
case  of  space-based  data,  also  the  full  range  of  latitudinal  phase  angles.  Assuming  a  new  Brightness  Data  point  was 
collected  every  90  minutes,  there  are  '-3 000  data  points  collected  during  the  nightly  operations  each  year  that  sample 
the  entire  accessible  range  space  for  the  vector  directions  of  solar  illumination  and  sensor  observations  for  a  satellite. 
The  record  of  Brightness  Data  represents  SSA  information. .  .hidden  in  plain  sight. 

References  1-3  utilize  the  two-facet  model  to  formulate  the  governing  equations  in  terms  of  two  satellite -attached 
reference  frames  (Section  2).  First  is  attached  to  the  satellite  body  (Body)  and  the  second  to  its  solar  panels  (Panels). 
The  Body  points  to  nadir  and  the  Panels  track  the  sun.  The  vector  directions  for  solar  illumination  and  sensor 
observations  are  computed  separately  for  the  Body  and  the  Panels.  As  a  result,  the  Brightness  Data  is  interpreted  in 
terms  of  entities  intrinsic  to  the  bidirectional  reflectance  distribution  functions  (BRDF)  for  the  Body  and  the  Panels. 
Furthermore,  the  use  of  the  principle  of  material  frame  indifference  allows  a  common  procedure  for  the  inversion  of 
data  collected  by  ground-base  and  space-based  sensors.  The  inversion  solves  for  pose-dependent  albedo-area  values 
for  the  satellite.  Body  and  Panels,  which  can  be  reused  in  order  to  predict  the  expected  value  of  Brightness  for  future 
observations.  Or  in  summary,  the  two  facet  model  is  used  in  the  Inversion  Model  as  well  as  the  Predictive  Model. 
This  is  fundamental  for  the  statistical  assessment  as  described  in  the  following  sections. 
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Figure  1.1:  An  idealized  schematic  of  synoptic  search.  The  West  and  East  designations 
for  the  Lambertian  regions  refer  to  the  direction  of  illumination. 

The  second  reason  can  be  addressed  by  using  the  statistical  assessment  procedure  described  in  this  paper.  It  draws 
upon  the  Bayes  network  approach  described  in  Reference  4  where  each  value  of  a  new  Brightness  Data  point  is 
compared  with  historical  values  of  Brightness  Data  points  under  matching  conditions  of  observation  geometry  in 
order  to  assess  if  the  satellite  Brightness  is  nominal  (NOM)  or  anomalous  (ANOM).  A  satellite  is  NOM  from  the 
viewpoint  of  photometry  if  its  observed  Brightness  matches  with  its  expected  Brightness  based  on  the  historical 
data.  A  satellite  is  ANOM  otherwise.  Such  assessment  is  performed  on  an  ongoing  basis  in  order  to  visualize 
continuity  in  the  propagated  Bayesian  belief,  where  belief  is  defined  as  marginal  probability  [Reference  5]. 

The  use  of  historical  data  for  assessment  of  NOM  and  ANOM  conditions  for  a  satellite  presents  a  peculiar  challenge. 
For  example,  consider  that  a  satellite  Brightness  may  be  denoted  as  ANOM  due  to  reasons  such  as  cross-tagging 
with  another  satellite,  visual  conjunction,  change  in  the  offset  angle  of  its  solar  panels,  180°  turn  (yaw  rotation)  in 
order  to  dump  momentum,  becoming  unstable,  appearance  of  unexpected  objects  in  the  sensor  field-of-view,  etc. 
[References  7-12].  However,  for  example,  cross-tagging  arises  due  to  incorrect  association  and  the  satellite  itself  has 
not  become  anomalous.  Similarly,  a  change  in  offset  angle  or  the  180°  turn  are  a  part  of  normal  operations  and  not 
anomalous  behavior.  Once  the  cross-tag  is  corrected,  the  new  offset-angle  is  computed,  or  the  180°  turn  is 
recognized,  the  satellite  status  would  revert  back  to  NOM,  although  it  would  be  new  NOM  for  which  there  may  or 


AMOS  2014  Technical  Conference 


may  not  be  sufficient  historical  data  available.  Also  during  the  period  over  which  the  historical  data  was  gathered,  , 
if  the  satellite  was  cross-tagged,  changed  its  solar  panel  offset,  etc.,  the  new  baseline  would  become  inaccurate. 

In  order  to  allow  the  changing  of  a  satellite  status  from  NOM  to  ANOM  and  vice  versa,  it  is  essential  to  possess  a 
consistent  baseline  which  may  be  used  to  compute  the  expected  value  of  Brightness  for  new  observations.  Such  a 
baseline  can  be  provided  by  the  Photometry  Model.  This  model  can  perform  inversion  of  existing  Brightness  Data 
and  then  use  the  computed  albedo-area  characteristics  to  predict  the  satellite  Brightness  for  new  observations  that 
are  collected  for  the  directions  of  illumination  and  sensor  line  of  sight  that  may  differ  from  those  in  the  existing  data. 

Marginal  probability  or  belief  is  the  probability  that  an  outcome  for  a  random  variable  equals  a  specific  value 
irrespective  of  the  outcome  of  the  other,  correlated  random  variables  IReference  6].  For  example,  for  jointly 
distributed  random  variables  X  and  F,  the  marginal  probability  P{X  =  x)  is  given  by: 

P{X  =  x)=  T.yP(X  =  X,  Y  =  y)  =  Y.yP(X  =  x\Y  =  y)P(Y  =  y) 

Or,  P(X  =  x)=  [P(X  =  x\Y)]  Equation  (1) 

Where  Ey  denotes  the  expected  value  or  expectation  (i.e.  mean)  with  respect  to  the  probability  distribution  of  Y.  Or, 
the  marginal  probability  for  P(X  =  x)  is  the  expected  value  for  the  probability  of  X  =  x  given  that  Y  =  y,  where  all 
possible  values  of  Y  are  covered  and  weighted  according  to  their  probabilities  and  ymin  ^  Y  ^  Ymax- the  context 
of  denoting  the  satellite  status  as  NOM  or  ANOM,  it  is  intrinsic  to  the  satellite  itself  and  is  independent  of  its 
directions  of  illumination  and  the  sensor  line  of  sight  (i.e.,  its  observation  geometry)  during  the  collection  of 
Brightness  Data.  Thus,  the  belief  that  a  satellite  is  NOM  or  ANOM  is  the  mean  probability  of  NOM  or  ANOM 
given  any  observation  geometry.  However,  there  are  infinite  combinations  of  observation  geometry  and  the 
Brightness  Data  collection  is  sparse.  Thus  a  compromise  is  to  compute  the  evolving  marginal  probability  (i.e.  belief 
propagation)  for  the  satellite  status  in  near-real  time,  as  each  new  point  of  Brightness  Data  becomes  available. 

2.  CHARACTER  OF  BRIGHTNESS  DATA 

Figure  2a  illustrates  a  notional  cluster  of  geosynchronous  satellites  (GEO  cluster)  and  how  its  size  compares  with  a 
notional  field  of  view  (FOV)  for  a  ground  or  space-based  sensor  that  is  utilized  for  synoptic  search.  A  GEO  cluster 
is  defined  as  a  collection  of  closely-spaced  satellites  that  are  station  kept  within  1°  angular  subtense  of  the  sky  with 
respect  to  a  ground-based  observer  IReference  7].  Thus  the  size  of  a  typical  GEO  cluster  is  smaller  than  the  FOV  for 
a  synoptic  search  sensor  and  the  Brightness  Data  for  the  entire  GEO  cluster  can  be  obtained  at  the  same  time. 

If  a  space-based  sensor  in  LEO  is  to  visit  each  GEO  cluster  once  in  each  pass,  the  mean  spacing  between  the 
Brightness  Data  points  for  a  GEO  cluster  will  be  roughly  90  minutes.  In  Figure  2.1b,  the  orbital  passes  are  denoted 
with  symbol  k  and  the  temporal  spacing  between  pass  (k-1)  and  k  is  denoted  as  Atk.  The  values  for  the  nightly  time 
spacing  Atk,  Atk+i  etc.  are  expected  to  be  unequal  due  to  differences  in  sensor  tasking  between  passes.  There  is  a  long 
daytime  gap  in  Brightness  Data  that  spans  between  the  morning  and  evening  terminators.  Figure  2.2  illustrates 
nightly  phase  angle  sampling  over  three  days.  For  a  ground-based  sensor,  such  Brightness  Data  could  be  collected 
with  a  back-and-forth  sweeping  motion  of  the  telescope  whereby  the  temporal  spacing  between  successive 
Brightness  Data  points  varies  as  a  function  of  the  relative  position  of  a  GEO  satellite  with  respect  to  the  sensor. 

Among  the  various  satellites  in  a  cluster,  denoted  as  Cj,  one  or  more  satellites  may  utilize  a  bus  type  different  from 
others  and  would  present  themselves  with  Brightness  values  that  differ  significantly  from  other  satellites  IReference 
7-8].  In  this  regard,  a  difference  in  excess  of  0.3  visual  magnitudes  or  2%  of  nominal  Brightness  may  be  considered 
significant  IReferences  10-13].  It  is  easier  to  obtain  correct  tagging  for  such  satellites.  References  7-8  discuss 
potential,  human  inspection  method  for  differentiating  between  such  satellites  with  as  few  as  3  to  5  points  of  data  as 
long  as  the  three  salient  regimes  of  Brightness  data  are  sampled;  namely  the  west  Lambertian  region,  the  central 
specular  region,  and  the  east  Lambertian  region.  The  west  Lambertian  region  occurs  in  the  hours  that  follow 
immediately  after  the  evening  terminator.  The  east  Lambertian  region  occurs  in  the  hours  prior  to  the  morning 
terminator.  The  Brightness  in  these  regions  is  dominated  by  diffuse  reflection  from  the  Body  and  the  Panel.  The 
central  specular  region  is  closer  to  midnight  and  the  signature  is  dominated  by  specular  reflection  from  the  Panel.  If 
a  nightly  collection  was  performed  once  every  90-minutes,  approximately  eight  new  data  points  would  be  collected 
per  night,  with  two  to  three  points  each  in  the  specular  region  and  in  the  west  and  east  Lambertian  regions. 
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Figure  2.1:  [a]  Relative  size  of  a  notional  GEO  cluster  and  sensor  field  of  view  (FOV) 
[b]  Notional  collection  of  one  set  of  Brightness  Data  per  orbital  pass  or  sweep 
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Figure  2.2:  Illustration  of  phase  angle  sampling  over  three  consecutive  days 

If  the  Brightness  and/or  proximity  for  satellites  Ci  and  Cj  are  close,  these  satellites  are  denoted  as  cluster  peers 
[Reference  12].  There  is  greater  potential  for  cross-tagging  between  such  cluster  peers.  There  can  also  be  a  flip-flop 
in  the  tagging  of  cluster  peers.  In  order  to  mitigate  such  cross -tagging,  the  availability  of  a  model  for  the  inversion 
and  prediction  of  Brightness  Data  is  important  in  order  to  support  the  calculations  for  statistical  assessment. 

3.  INVERSION  AND  PREDICTION  OF  BRIGHTNESS  DATA 

Due  to  the  variable  sensor  tasking,  orbital  motion  of  the  space-based  sensor  and  a  continuous  change  in  the  solar 
declination  angle,  the  observation  geometry  is  different  for  each  pass.  Also,  the  satellite  operators  periodically 
change  the  solar  panel  offset  angle  and  may  also  subject  the  spacecraft  to  180°  turn  in  order  dump  momentum. 
Furthermore,  the  aging  of  satellite  materials  reduces  the  Brightness  of  the  Body  and  Panel.  As  a  result,  it  is 
necessary  to  maintain  an  updated  baseline  for  the  nominal  Brightness  for  the  satellite.  Such  a  baseline  can  be 
extracted  by  the  inversion  of  Brightness  Data  using  the  Photometry  Model  reported  in  References  1-3.  In  this  model, 
the  satellite  is  represented  with  a  two  facet  model  and  the  data  inversion  is  posed  to  solve  for  the  intrinsic  albedo - 
areas  for  the  Body  and  the  Panel,  as  well  as  the  Panel  offset  angle.  The  term  ‘intrinsic’  is  meant  to  describe  the 
scaling  factors  that  appear  in  a  BRDF  relationship.  For  example,  Lambert’s  cosine  law  defines  the  change  of 


AMOS  2014  Technical  Conference 


Brightness  as  a  function  of  the  phase  angle.  The  numerical  value  of  Brightness  contribution  by  the  Body  and  Panel 
under  Lambertian  conditions  is  given  by  the  multiplication  of  the  cosine  term  with  their  respective  albedo-area 
values.  Once  the  albedo-area  and  offset  angle  are  known,  it  is  feasible  to  predict  the  satellite  Brightness  for  the 
observation  conditions  in  the  future  passes  using  the  analytical  expressions  of  the  BRDF  for  the  Body  and  Panel. 


In  order  to  provide  a  continuously  updated  baseline,  the  Photometry  Model  application  is  performed  using  a  time 
slider  approach  illustrated  in  Figure  3.1a.  The  origin  of  the  time  slider  is  positioned  at  pass  k  =  0,  where  it  is  deemed 
prior  knowledge  that  the  satellite  is  NOM  for  sufficient  number  of  prior  passes.  In  terms  of  the  notation  of  Figure 
3.1a,  the  orbital  passes  that  comprise  prior  data  are  assigned  values  of  k  less  than  zero.  The  new  data  is  when  k  is 
greater  than  or  equal  to  zero.  The  prior  data  is  utilized  for  the  statistical  characterization  of  NOM  behavior  of  the 
satellite.  To  this  end,  it  is  sufficient  if  the  number  of  prior  passes  provide  a  minimum  of  thirty  points  of  Brightness 
Data  for  each  of  the  three  salient  regions;  namely  the  west  and  east  Lambertian  region  and  the  specular  region.  At  a 
median  rate  of  two  to  three  points  of  data  per  day  for  the  specular  and  two  Lambertian  regions,  the  number  of  passes 
in  two  weeks  is  expected  to  provide  sufficient  data  to  perform  inversion  using  the  two -facet  model.  The  Photometry 
Model  is  used  for  the  inversion  of  the  prior  data  (k  <  0)  and  for  the  prediction  of  Brightness  for  the  new  data  (k  >  0). 


[a] 


Time  Slider,  k  =  0 


NOM 


Prior  Data  k  <  0  k  >  0  New  Data 


•••  ^ 


•  •• 


Model  is  used  for  data  inversion 


Model  is  used  for  prediction 


Figure  3.1:  [a]  Model  solution  from  data  inversion  on  prior  data  is  used  to  predict  brightness  for  new  data 

[b]  GEO  satellite  body  points  to  nadir  and  solar  panels  track  the  sun.  Offset  angle  is  nonzero 

[c]  Model  extracts  Body  and  Panel  albedo  areas  and  offset  angle  for  the  Panel  from  prior  data 


In  order  to  determine  the  location  of  the  time  slider,  it  is  useful  to  consider  that  the  statistical  assessment  is  an 
ongoing  process.  In  order  to  initialize  the  assessment,  an  analyst  needs  to  mark  a  Brightness  Data  point  as  pass  k  =  0 
such  that  the  satellite  is  NOM  prior  to  this  pass  for  a  minimum  of  two  weeks.  Once  the  statistical  assessment  begins, 
the  time  slider  is  moved  forward  to  a  temporal  position  at  which  the  level  of  belief  for  NOM  status  for  a  satellite 
fulfills  a  threshold  limit  for  the  level  of  confidence  [Reference  23].  This  procedure  is  described  further  in  Section  15. 
However,  this  will  result  in  a  different  location  for  the  time  slider  for  each  satellite  in  the  cluster.  Thus,  it  is  simpler 
to  move  the  time  slider  only  to  the  largest  pass  number  for  which  all  satellites  in  the  cluster  are  NOM. 
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Note  that  just  as  the  Brightness  Data  differs  from  one  satellite  to  another,  the  Photometry  Model  accuracy  for 
inversion  and  prediction  also  varies  from  one  satellite  to  another  and  from  one  type  of  observation  geometry  to 
another.  Indeed,  unless  proven  otherwise,  it  is  useful  to  assume  that  the  Photometry  Model  is  imperfect  and  has 
limited  accuracy.  Thus,  a  dimensionless  ratio  Vjk  is  utilized  in  order  to  quantify  the  difference  between  the  observed 
and  predicted  Brightness  values.  This  ratio  is  denoted  as  Brightness  Ratio.  Specifically: 

=  Vj^  —  1  Equation  (2a) 

^ik  ~  Equation  (2b) 

^mjk 

where  is  the  Brightness  observed  for  satellite  j  and  pass  k,  and  is  the  Brightness  predicted  by  the  model  for 
satellite  j  and  pass  k.  Thus,  if  the  model  was  a  good  predictor  for  satellite  j,  rjk  would  be  a  zero  mean  random 
variable.  A  probability  mass  function  for  rjk  can  now  be  constructed  using  prior  data  (i.e.  k  <  0).  Since  is 

continuous,  a  continuous,  real  random  variable  R  is  defined  such  that  its  range  spans  all  feasible  values  of 

4.  ASSUMPTIONS,  INPUT  DATA  REQUIREMENTS  AND  OUTPUT 

The  following  assumptions  are  utilized  in  this  work: 

1)  The  NOM  or  ANOM  state  of  each  satellite  is  independent  of  other  satellites,  except  for  cluster  peers. 

2)  Cross-tagging  is  feasible  only  between  cluster  peers. 

3)  The  probability  that  a  satellite  is  NOM  or  ANOM  during  pass  (k-Fl)  depends  only  on  its  state  during  pass  k 
and  not  on  its  state  during  the  passes  prior  to  k.  The  probability  that  a  satellite  is  NOM  or  ANOM  during 
pass  k  depends  only  on  its  state  during  pass  (k-1)  and  not  on  its  state  during  the  passes  prior  to  (k-1),  etc. 

4)  Each  point  of  observation  data  for  all  satellites  in  a  cluster  is  taken  at  the  same  time 

Eour  sets  of  input  data  are  required  as  follows: 

1)  The  observations  data  is  available  for  all  satellites  in  a  cluster  as  soon  as  it  becomes  available  (i.e.,  in  near 
real  time). 

2)  The  following  metadata  is  provided  along  with  each  point  observation  data:  Julian  date.  Brightness,  the 
coordinates  for  the  satellite,  sensor  and  the  sun. 

3)  The  long  run  probability  that  a  satellite  is  NOM  at  any  time,  Long  run  probability  means  long  term 
probability  or  probability  of  NOM  irrespective  of  time,  which  is  the  same  as  marginal  probability  or  belief. 
This  value  is  estimated  by  the  fraction  of  passes  in  the  historical  data  that  spans  an  extended  duration  (long 
run  historical  data)  where  the  state  of  the  satellite  is  identified  as  NOM. 

4)  The  terminator-to-terminator  transition  probability  for  a  satellite.  This  is  probability  that  the  satellite  is 

NOM  at  the  evening  terminator  given  that  satellite  was  NOM  at  the  morning  terminator  (e.g.,  this 
probability  is  lower  for  cluster  peers  that  are  prone  to  cross-tag).  It  is  denoted  by  where  the  subscript 

Gap  is  in  reference  to  the  lack  of  data  during  this  extended  period  of  time.  This  is  estimated  by  the  fraction 
of  days  in  the  long  run  historical  data  when  the  satellite  is  denoted  as  NOM  at  the  last  pass  for  the  night  at 
the  morning  terminator  and  then  denoted  as  NOM  at  the  first  pass  for  the  following  evening  terminator. 

Note  that  the  estimated  values  for  long  run  probability  and  terminator-to-terminator  transition  probability 
are  based  on  prior  or  historical  data.  They  are  defined  by  the  user.  They  can  be  season  or  age  dependent  for  each 
satellite.  They  can  be  reset  by  the  user  in  the  midst  of  a  long-term,  continuous  statistical  assessment  whereby  the 
updated  values  are  utilized  for  all  points  of  Brightness  Data  thereafter. 

The  output  comprises  a  numerical  value  for  a  measure  of  belief  in  the  NOM  or  ANOM  state  of  each  satellite  in  a 
GEO  cluster.  This  output  is  expected  to  become  available  within  minutes  after  the  receipt  of  each  new  point  of 
Brightness  Data.  The  utility  of  the  ANOM  notification  is  to  inform  a  user  that  the  statistical  assessment  has  detected 
evidence  that  the  Brightness  Ratio  is  different  from  when  the  satellite  was  NOM.  The  ANOM  notification  can  also 
include  a  likely  cause  for  this  difference. 


AMOS  2014  Technical  Conference 


5.  PRIOR  PROBABILITY  DISTRIBUTION  FUNCTIONS  FOR  THE  BRIGHTNESS  RATIO 


The  probability  distribution  function  (PDF)  for  the  Brightness  Ratio,  calculated  using  the  prior  data  (or,  prior  PDF) 
is  needed  as  baseline  for  the  assessment  of  occurrence  of  ANOM  status  in  the  new  data  (Figure  3.1).  This 
calculation  is  at  the  heart  of  the  statistical  assessment.  The  prior  PDF  is  required  for  both  NOM  and  ANOM 
conditions.  A  procedure  for  the  calculation  of  the  NOM  and  ANOM  PDFs  is  described  in  the  following: 

Since  the  satellite  is  deemed  NOM  prior  to  the  placement  of  the  time  slider,  the  prior  PDF  for  NOM  can  be 
calculated  using  the  Brightness  Ratio  data  for  passes  k  <  0.  A  continuous  two  week  period  for  the  prior  data  is 
expected  to  provide  over  a  hundred  values  of  NOM  data  for  Brightness  Ratio.  These  values  will  be  distributed  over 
the  west  and  east  Lambertian  and  specular  regions.  A  longer  period  of  prior  data  could  be  considered,  but  it  is  not  a 
requirement.  This  is  because  there  are  limits  on  how  long  a  continuous  stretch  of  NOM  prior  data  can  be  feasible  on 
a  routine  basis.  For  example,  if  two  cluster  peer  satellites  can  become  cross-tagged  frequently,  the  mean  duration  of 
time  when  these  satellites  are  NOM  is  shorter.  As  a  result  of  limiting  to  a  relatively  short  prior  region,  the  prior  PDF 
for  the  NOM  state  is,  in  reality,  closer  to  a  probability  mass  function  (PMF),  which  is  estimated  by  constructing  a 
histogram  for  the  values  of  the  Brightness  Ratio.  Also, 

/(  ^;fc)  =  /(  rjk  I  NOMj  )  P{  NOMj)  +  /(  rj„  \  ANOMj)  P(ANOMj')  Equation  (3) 


Where 

/( 7}/^)  is  the  PDF  for  the  Brightness  Ratio  for  pass  k  for  satellite  Q 

/( rjk  I  NOMj  )  is  the  PDF  for  the  Brightness  Ratio  for  pass  k  given  the  satellite  Q  is  NOM, 

/( rjk  I  ANOMj)  is  the  PDF  for  the  Brightness  Ratio  for  pass  k  given  the  satellite  Q  is  ANOM, 

P(  NOMj)  is  the  long  run  probability  that  the  satellite  is  NOM,  also  denoted  as  tt^,  and 
P(ANOMj)  is  the  long  run  probability  that  the  satellite  is  ANOM,  also  denoted  as 
pIaNOMj)  =  1  -  tTjv 

The  calculation  of  PDF  for  rjk  given  ANOM  is  more  detailed  because  ANOM  is  defined  relative  to  current  state  of 
NOM.  Also,  ANOM  is  a  collective  term  for  a  variety  of  conditions.  Assuming  that  these  conditions  are  mutually 
disjoint, 

P(^ANOMj)  =  P{p f f  setChange j)  +  P{CrossTagj)  +  P(  Turrij)  +  P(  Unstablej)  +  P(  UCTj)  ... 


Equation  (4a) 


and  P(  OffsetChangej  \  ANOMj) 


P{OffsetChangej) 
P  (ANOMj) 


,  P(CrossTagj  \  ANOMj)  = 


P(CrossTag  j) 
P  (ANOMj) 


,  etc.  Then, 


f{rjj^  I  ANOMj)  =  f{rjj^  \  OffsetChangej)  P(OffsetChangej  \  ANOMj) 
+  f{rjj^  I  CrossTagj)  P{CrossTagj  \  ANOMj) 

+ /( ^jk  I  TurUj  )  P(Turnj  \  ANOMj) 

+  /(  Tjj^  I  Unstablej )  P(Unstablej  \  ANOMj) 

+  /(  Tjj,  I  UCTj)  P{UCTj  I  ANOMj)  . . . 


Or, 


Equation  (5  a) 

f^Vjj^  n  ANOMj)  =  f^Tjj^  I  0 ffsetChangej)  P (^OffsetChangej) 

+  /(^7c  I  CrossTagj)  P(CrossTagj)  +  f{rjj^  \  Turnj)  P(Turnj) 

+  /(^77c  I  Unstablej)  P(Unstablej)  +  f{rjj^  \  UCTj)  P(UCTj)  ... 

Equation  (5b) 
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Where, 


P{OffsetChangej)  is  the  long  run  probability  of  change  in  the  panel  offset  angle  for  satellite  Q 
P{Crosstagj)  is  the  long  run  probability  of  cross-tag  for  Q 
P(  Turrij)  is  the  long  run  probability  of  a  180°  turn  for  Q 
P(  Uns table j)  is  the  long  run  probability  that  satellite  Q  is  unstable 

P(  UCTj)  is  the  long  run  probability  of  a  mistag  for  Q  due  to  the  appearance  of  UCT  in  the  sensor  FOV 

And, 

/( rjk  I  OffsetChangCj)  is  the  PDF  of  Brightness  Ratio  when  the  Q  panel  offset  is  modified 
/( rjk  I  CrossTagj)  is  the  PDF  of  Brightness  Ratio  when  Q  is  cross-tagged  with  a  peer 
/( rjk  I  Turrij)  is  the  PDF  of  Brightness  Ratio  when  Q  is  turned  180° 

/( rjk  I  Unstablej)  is  the  PDF  of  Brightness  Ratio  when  Cj  is  unstable 

/( rjk  I  UCTj)  is  the  PDF  of  Brightness  Ratio  when  a  mistag  is  caused  by  an  unexpected  UCT 

The  long  run  probability  values  for  the  anomalous  conditions  can  be  identified  from  the  Brightness  Ratios  for  the 

historical  data  because  each  condition  has  a  tell-tale  feature.  In  order  to  locate  these  features,  it  is  necessary  to 

consider  plots  of  Brightness  Ratio  versus  phase  angle  for  each  night  of  historical  data.  These  plots  can  be  generated 

automatically  and  visualized  in  order  to  locate  specific  anomalies  based  on  their  features  as  described  below: 

•  Change  in  the  solar  panel  offset  causes  the  Brightness  Ratio  to  become  anomalous  only  in  the  specular  region. 
The  Brightness  Ratio  in  the  west  and  east  Lambertian  regions  is  minimally  affected  because  the  Panels  are 
highly  specular  and  the  effect  of  offset  angle  change  fades  significantly  at  high  phase  angles.  Thus 
P{OffsetChangej)  equals  the  fraction  of  passes  in  the  historical  data  where  the  Brightness  Ratio  is  ANOM 
only  in  the  specular  region. 

•  Occurrence  of  cross-tag  causes  the  Brightness  Ratios  for  cluster  peers  to  become  anomalous  simultaneously. 
Thus,  P (^CrossTagj)  equals  the  fraction  of  passes  in  the  historical  data  where  the  Brightness  Ratio  has  become 
simultaneously  anomalous  for  the  cluster  peers.  It  would  also  be  feasible  to  verify  the  cross-tag.  For  example, 
let  satellite  Ci  and  Cj  be  the  cluster  peers  that  cross-tag.  Then  the  Brightness  Ratio  for  Ci  could  be  recomputed 
using  the  term  Imjk  instead  of  7^^/^  and  vice  versa  for  Q.  This  would  render  the  Brightness  Ratio  NOM  if  the 
peers  were  cross-tagged. 

•  The  effect  of  180°  yaw  turn  by  a  satellite  causes  the  sides  of  the  Body  to  switch  with  respect  to  the  illumination 
by  sunlight  as  illustrated  in  Figure  5.1.  Imagine  that  the  satellite  is  shown  in  a  pass  prior  to  the  180°  yaw  turn. 
The  normal  to  one  of  the  body  facets  points  eastwards.  Upon  the  180°  yaw  turn,  this  normal  will  point 
westwards.  Thus,  the  ANOM  status  resulting  from  a  180°  yaw  turn  causes  the  values  of  Brightness  Ratio  to 
switch.  The  values  that  normally  occurred  in  the  west  Lambertian  region  will  now  appear  in  the  east  Lambertian 
region,  and  vice  versa  (Figure  1.1).  Thus  P(  Turrij)  equals  the  fraction  of  passes  in  the  historical  data  where  the 
Brightness  Ratio  is  ANOM  in  a  manner  that  depicts  an  east-west  switch. 

•  Effect  of  gross  instability  is  loss  of  attitude  control.  This  creates  opportunities  for  the  occurrence  of  specular 
glints  in  the  west  and  east  Lambertian  regions  and  the  occurrence  of  diffuse  reflection  conditions  in  the  specular 
regions.  Thus,  the  ANOM  status  for  gross  instability  is  plausible  if  the  values  of  high  Brightness  are 
uncorrelated  to  the  phase  angle.  Thus  P(  Unstablej)  equals  the  fraction  of  passes  in  the  historical  data  where 
the  Brightness  Ratio  is  ANOM  in  a  manner  that  depicts  lack  of  correlation  with  the  phase  angle. 

•  A  orbital  debris  in  LEO  or  a  higher  orbit  with  its  velocity  vector  aligned  with  the  view  direction  for  a  ground  or 
space-based  sensor  can  cause  presence  of  a  UCT  (i.e.  an  extra  object  in  close  proximity  of  a  GEO  cluster), 
potentially  resulting  in  a  mistag  when  the  object  association  is  performed  based  only  on  metric  data.  The  same 
UCT  is  unlikely  to  be  present  the  next  time  the  sensor  visits  the  GEO  cluster  and  the  error  due  to  incorrect 
association  may  or  may  not  be  rectified.  Indeed,  another  UCT  may  appear  subsequently  and  propagate  the  cycle 
of  incorrect  association  further  [Reference  10].  Thus  a  situation  where  the  Brightness  Ratio  is  anomalous  in  one 
pass  but  is  nominal  in  the  subsequent  passes  may  be  due  to  a  UCT  that  appears  momentarily  in  the  sensor  field 
of  view.  Thus  P(  UCTj)  equals  the  fraction  of  passes  in  the  historical  data  where  the  Brightness  Ratio  is  ANOM 
in  a  manner  that  depicts  an  isolated  anomaly. 

The  calculation  procedure  for  the  conditional  PDFs  (e.g.  /(  rjj^  \  0 f  f  setChange j))  is  described  in  Section  7. 
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Orbit  Normal 


Figure  5.1:  Effect  of  an  180°  yaw  turn  on  the  orientation  of  a  GEO  satellite 

6.  BAYES  NETWORK 


Figure  6.1:  Flow  of  information  within  the  Bayes  network  for  satellite  Cj 

Figure  6.1  shows  a  graphical  model  for  the  Bayes  network  to  perform  statistical  assessment  of  NOM  and  ANOM 
status  for  satellite  Q  using  Brightness  Data.  This  graphical  model  builds  upon  the  Bayes  network  reported  in 
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Reference  4.  As  in  Reference  4,  it  contains  GEO  status  notes  (blue  nodes)  that  receive  Brightness  evidence  (orange 
nodes)  after  each  pass.  The  arrows  in  the  network  have  a  direction  but  they  do  not  loop  around.  In  other  words,  it  is 
a  directed  acyclic  graph.  The  network  is  directional  in  that  it  defines  the  qualitative  dependencies  between  the  nodes 
[Reference  5].  The  number  of  GEO  status  nodes  is  equal  to  the  number  of  passes  performed  subsequent  to  the 
placement  of  time  slider  at  k  =  0.  When  the  time  slider  is  moved  forward  to  a  new  position,  the  Bayes  Network  is 
reset  (Sections  15-16).  The  GEO  status  node  after  pass  k  is  the  parent  node  for  the  GEO  status  after  pass  k-^l  and  the 
dependency  between  nodes  is  defined  in  terms  of  conditional  probability.  The  temporal  spacing  between  the  nodes 
is  non-uniform,  which  corresponds  to  the  time  elapsed  between  the  successive  collection  of  Brightness  Data  points. 

In  Eigure  6.1,  P  (  NOMjj^  \  )  is  the  conditional  probability  that  the  status  of  satellite  Cj  is  NOM  at  pass  k  given 

the  Brightness  Ratio  for  the  data  collected  in  pass  k.  The  calculation  of  this  conditional  probability  is  performed 
using  Bayes  theorem.  This  calculation  is  described  in  Section  7.  The  quantity  P  (  NOMjj^  \  Vjj^  )  occurs  often  in  the 
following  calculations  and  it  is  given  a  simpler  notation,  Xjk.  Note  that  the  notation  rjj^  is  meant  to  imply  that  the 
value  of  the  continuous  random  variable  R  resides  in  an  infinitesimal  neighborhood  of  a  given  number  rjj^.  This  is 
because  R  is  a  continuous  random  variable  and  thus  the  probability  that  R  exactly  equals  any  given  number  Tjj^  is 
zero.  This  consideration  is  absent  in  the  case  of  the  probability  of  NOM  or  ANOM  states  because  they  are  discrete. 

In  Eigure  6.1,  P  |  is  the  conditional  probability  that  the  status  of  satellite  Cj  is  NOM  at  the  time 

of  collection  of  data  in  pass  (k-^l)  given  satellite  Cj  status  is  NOM  at  the  time  of  collection  of  data  in  pass  k.  The 
term  P  |  ANOMjj^)  is  the  conditional  probability  that  the  status  of  satellite  Cj  is  NOM  at  the  time  of 

collection  of  data  in  pass  (k-Fl)  given  satellite  Cj  status  is  ANOM  at  the  time  of  collection  of  data  in  pass  k.  These 
two  terms  define  the  probability  of  transition  of  the  NOM  state  during  the  period  when  no  new  data  is  collected. 
These  two  terms  are  the  first  column  of  the  transition  probability  matrix,  which  is  computed  using  the  Markov  chain. 
This  calculation  is  described  in  Section  8.  It  uses  a  simpler  notation  for  the  conditional  probability  terms  as  follows: 

N^'k+i  —  P  I  Equation  6a 

=  p  I  ANOMjk)  Equation  6b 

This  notation  uses  two  subscripts  and  two  superscripts.  The  left  and  right  superscripts  denote  the  state  and  pass 
number,  k.  The  left  and  right  subscripts  denote  the  state  and  pass  number,  Two  additional,  related  terms  are 
defined  in  Equation  6  and  they  appear  in  relation  to  the  Markov  matrix  calculations  in  Section  8. 

aTu+i  =  1  -  nTu+i  Equation  6c 

=  1  -  Equation  6d 

The  flow  of  information  in  the  network  can  be  visualized  using  the  GEO  status  node  for  pass  (k-i-1).  The  transition 
probability  terms  are  utilized  to  compute  the  probability  that  the  satellite  is  NOM  prior  to  the  collection  of  data  in 
pass  (k-Fl)  based  on  the  probability  that  the  satellite  is  NOM  for  pass  k.  The  evidence  of  Brightness  Ratio  for  the 
data  collected  in  pass  (k-^l)  is  utilized  to  compute  the  probability  that  the  satellite  is  NOM  after  the  data  has  been 
collected  in  pass  (k-^l).  Then  the  transition  probability  terms  are  utilized  to  compute  the  probability  that  the  satellite 
is  NOM  prior  to  the  collection  of  data  in  pass  (k-\-2);  etc.  This  calculation  is  performed  using  the  Chapman- 
Kolmogorov  equation  and  the  Markov  chain,  and  it  is  described  in  Section  9. 

In  following  the  Bayes  Network  directed  graph,  the  numerical  value  of  probability  that  satellite  Cj  is  NOM  evolves 
with  the  collection  of  each  new  point  of  Brightness  Data.  The  value  of  Brightness  in  each  pass  is  a  function  of  the 
observation  geometry,  namely  the  directions  of  illumination  and  the  sensor  line  of  sight.  Thus,  the  conditional 
probability  that  satellite  Cj  is  NOM  after  pass  k  given  the  Brightness  Ratio  rjk  can  also  be  interpreted  as  the 
conditional  probability  that  satellite  Cj  is  NOM  after  pass  k  given  the  observation  geometry  for  pass  k. 

The  probability  density  function  (PDE)  for  the  Brightness  Ratio  can  be  quite  different  from  one  satellite  to  another, 
even  if  the  satellites  are  cluster  peers.  This  is  due  to  limitation  on  accuracy  that  can  be  achieved  by  the  Inversion 
Model.  An  empirical  description  of  this  PDE  is  computed  using  the  data  for  passes  k  <  0  and  it  is  used  subsequently 
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to  calculate  Xjk  for  passes  k  >  0.  In  the  case  of  a  cross-tag  between  cluster  peers,  the  values  of  Xjk  get  affected 
simultaneously  for  both  peer  satellites.  In  order  to  detect  such  simultaneous  changes,  it  is  beneficial  to  transform  the 
values  of  Xjk  to  z-score  in  order  to  support  the  calculation  of  a  measure  of  belief  and  for  hypothesis  testing  to  assess 
the  probabilities  of  Type  I  and  Type  II  errors  [References  23  and  26]. 

In  order  to  assign  a  numerical  value  to  the  belief  placed  in  the  status  of  satellite  Cj,  we  need  know  the  conditional 
probability  irrespective  of  the  observation  geometry  (see  Equation  1),  which  is  the  marginal  probability  that  satellite 
Cj  is  NOM.  In  order  to  compute  such  marginal  probability,  we  need  to  sample  the  entire  range  space  of  observation 
geometry  at  a  sufficient  number  of  locations.  This  is  mostly  not  feasible  because  the  observation  geometry  is 
governed  by  automated  tasking  for  synoptic  search  by  the  sensor.  As  a  result,  it  is  necessary  to  perform  an  iterative 
assessment  of  belief  based  on  a  limited  number  of  data  points  for  each  satellite  in  a  cluster.  This  is  performed  using 
the  standard  z-score,  calculation  of  likelihood  rations  and  hypothesis  testing  (see  Sections  11  to  15). 

7.  BAYES  THEOREM 


The  purpose  of  the  calculations  described  in  this  section  is  to  determine  the  Brightness  evidence,  Xjk  (Figure  6.1).  In 
this  regard,  the  probability  that  satellite  Cj  is  NOM  after  pass  k  and  the  Brightness  Ratio  resides  in  an  infinitesimal 
neighborhood  of  rjj^  can  be  defined  using  the  expression  for  conditional  probability  in  two  ways  as  follows: 

P[  NOMjk  n  (rjk  <R  <  rjk  +  Arj^)  )  =  P(  NOMjk  \  (rj^.  <R<  rj^.  +  Arj^)  )  P(rjk  <R<  +  Arj^) 

Equation  7  a 

p[  NOMjy.  n  (r^fc  <R<  rjk  +  Avj^)  }  =  P{  (rj^  <R<  +  Arj^)  \  NOMj^}  P{  NOMjy.  ) 

Equation  7b 


Substituting  Equation  7a  into  Equation  7b  and  rearranging. 


P(  NOMjk  I  (rjk  <R<  rjk  +  Arj^)  ) 


P(  jrj,  <R<  rj,  +  Arj,)  \  NOMj,  )  P(  NOMj,  ) 
P(rjk  <R  <  rjk  +  Arjk) 


Equation  8 


Equation  8  is  now  re-written  using  the  probability  density  function  (PDF)  for  the  Brightness  Ratio  /(  ),  which  is 

suitably  interpreted  as  the  probability  of  R  taking  values  in  the  small  neighborhood  about  the  given  value: 


P(NOMj„  I  rj„)  = 


f{rji,  \N0Mji,)P(N0Mji,) 


fiTju)  =  /(  Vju  I  NOMj^  )  P(  NOMj^)  +  /(  \  ANOMj,  )  P(  ANOMj,  ) 


Equation  9a 
Equation  9b 


The  above  derivation  is  essentially  Bayes  theorem.  The  left  hand  side  of  Equation  9a  is  Xjk.  Equation  9b  expresses 
the  denominator  of  Equation  9a  using  the  law  of  total  probability.  The  calculation  of  P(  NOMj^  )  is  connected  to  the 
calculation  of  transition  probability  and  it  is  described  in  Sections  8  and  9. 

The  conditional  PDF  /(  |  NOMjj^  )  is  computed  as  empirical  PDF  from  the  prior  data,  which  is  denoted  by 

passes  k  <  0  (Figure  3.1).  Note  that  this  is  permissible  because  the  time  slider  is  positioned  such  that  the  satellite 
status  is  NOM  for  a  minimum  of  two  weeks  prior  to  the  location  k  =  0.  This  data  is  utilized  as  input  to  the  Inversion 
Model  and  for  the  calculation  of  Note  that  such  inversion  calculation  is  performed  upon  the  translation  of  the 
time  slider  to  each  new  position  so  that  the  empirical  PDF  /(  |  NOMjj^  )  may  evolve  with  the  passage  of  time. 


The  f(rjj^)  calculation  includes  contributions  from  NOM  and  ANOM  conditions  (Equation  3).  Note  that  the  first 
term  in  Equation  3  is  the  NOM  term,  which  is  the  same  as  the  numerator  of  Equation  9  for  pass  k  =  0.  The  second 
term  in  Equation  3  is  the  ANOM  term.  It  is  given  by  Equation  5a  by  considering  the  contributions  from  five  disjoint 
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anomalies.  Among  these  anomalous  states,  the  CrossTag,  Turn,  and  Unstable  are  Boolean  states.  They  are  either  true 
or  false.  The  OffsetChange  and  UCT  are  continuous  states.  The  changed  solar  panel  offset  may  be  any  real  number 
within  a  range  of  permissible  values.  The  Brightness  of  a  UCT  may  be  any  real  number,  although  it  will  typically  be 
a  fainter  object.  The  PDF  calculation  for  the  discrete  and  continuous  ANOM  states  is  described  in  the  following 
based  on  a  combination  of  arguments  derived  from  the  observations  data  and  physical  considerations. 

The  calculation  of  /(  Tjj^  \  CrossTag j)  is  performed  as  follows.  Consider  that  satellite  Cj  is  incorrectly  tagged  as  its 
cluster  peer,  denoted  as  Then  the  Brightness  Ratio  given  CrossTag  can  be  expressed  as: 

{  Imjk  /  ImTk  }  Vjk  Equation  10a 

^CrossTag  _  ^  Equation  10b 

In  Equation  10a,  the  value  of  given  by  Equation  2b.  It  is  scaled  by  the  ratio  of  predicted  Brightness  for  satellite 
Cj  to  the  predicted  Brightness  for  satellite  •  The  predicted  brightness  for  Cj  and  assume  both  satellites  to 
be  NOM.  The  resulting  values  of  are  used  to  compute  the  empirical  PDF  for  /(  Vjj^  \  CrossTag  j). 

The  calculation  of  /(  |  Turrij)  is  performed  as  follows.  Consider  that  satellite  Cj  has  undergone  a  180°  yaw  turn. 

The  turned  satellite  is  denoted  as  Then  the  Brightness  Ratio  given  Turn  can  be  expressed  as: 

^jurn  ^  !  furn  j  Equation  1  la 

rj^rn  ^  ^Turn  _  Equation  11b 

In  Equation  11a,  the  value  of  Vjj^  given  by  Equation  2b.  It  is  scaled  by  the  ratio  of  predicted  Brightness  for  satellite 
Cj  to  the  predicted  Brightness  for  satellite  The  rotation  of  the  satellite  modifies  the  Brightness  contributions 

of  the  Body  and  Panel  to  the  total  Brightness  for  a  satellite  as  described  using  Figure  5.1.  The  modified  contributions 
of  the  Body  and  Panel  can  be  computed  using  the  Predictive  Model.  The  modified  projected  view  upon  the  yaw  turn 
is  utilized  in  order  to  compute  the  predicted  brightness  The  resulting  values  of  r^^^^  are  used  to  compute  the 

empirical  PDF  for  /( 7}/^  |  Turrij). 

The  calculation  of  /(  rjj^  \  Unstablej)  is  performed  as  follows.  Consider  that  satellite  Cybecomes  unstable  and  it  is 
now  denoted  as  the  Brightness  Ratio  for  the  Unstable  condition  can  be  expressed  as: 

^unstable  =  |  |  Equation  12a 

^unstable  ^  ^Unstable  _  ^  Equation  12b 

In  Equation  12a,  the  value  of  Vjj^  given  by  Equation  2b.  Note  that  when  a  satellite  is  unstable,  its  projected  view  with 
respect  to  the  sensor  is  independent  of  the  solar  phase  angle.  Or,  its  Brightness  at  high  phase  angles  can  be  larger 
than  that  at  low  phase  angles.  Thus,  in  the  calculation  of  ^  the  value  of  is  randomly  distributed 

between  the  maximum  and  minimum  values  possible  for  the  Brightness  of  satellite  Cj  while  in  its  NOM  state.  These 
values  are  denoted  as  bounds,  which  are  computed  using  the  Predictive  Model.  The  range  for  the  resulting  values  for 
^Unstable  ^  function  of  the  phase  angle  because  is  a  single  number.  These  values  for  are  used  to 

compute  the  empirical  PDF  for  /(  |  Unstablej).  Note  that  this  PDF  also  becomes  a  function  of  the  phase  angle. 

The  calculation  of  /(  rjj^  \  OffsetChangej)  is  performed  as  follows.  Consider  that  satellite  exchanges  its  solar 
panel  offset  and  it  is  now  denoted  as  the  Brightness  Ratio  for  the  OffsetChange  condition 

can  be  expressed  as: 


CrossTag  _ 
jk 


CrossTag  _ 
'^jk 
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k;."""”'"  =  {  /„/,  /  )  %  Equation  13a 

=  nWW««.na  _  j  Equation  13b 

Note  that  when  a  satellite  changes  its  Panel  offset,  the  projected  view  of  the  Panel  with  respect  to  the  sensor  changes 
while  the  projected  view  of  the  body  remains  unchanged.  Since  the  Panels  are  highly  specular,  their  contribution  to 
the  satellite  Brightness  in  both  Lambertian  regions  is  small  as  compared  to  the  Body  (Figure  1.1).  Or,  the  value  of 
remains  largely  unchanged  in  the  Lambertian  regions  in  spite  of  the  change  in  the  offset  angle. 

In  the  specular  region,  randomly  distributed  between  the  maximum  and  minimum  values  of 

Brightness  that  can  occur  for  a  NOM  satellite  Cj  in  the  specular  region.  These  maximum  and  minimum  values  are 
computed  using  the  Predictive  Model.  The  range  for  the  resulting  values  for  ^  function  of  the 

phase  angle  in  the  specular  region  only.  These  values  for  compute  the  empirical  PDF  for 

/( rjk  I  OffsetChangej).  Note  that  this  PDF  is  a  function  of  the  phase  angle  during  pass  k  when  it  is  in  the 
specular  region.  The  PDF  for  /(  Tjj^  \  OffsetChangej)  for  the  Lambertian  regions  is  same  as  that  for  the  NOM 
state. 

The  calculation  of  |  UCTj)  is  performed  as  follows.  This  calculation  is  similar  to  the  Unstable  state  of 

ANOM.  Consider  that  a  UCT  is  incorrectly  tagged  as  satellite  Cj.  The  incorrectly  tagged  UCT  is  now  denoted  as 
Cj^cT ^  Then  the  Brightness  Ratio  for  the  UCT  condition  can  expressed  as: 

=  {  Imjk  /  lm%  }  ^jk  Equation  14a 

-  1  Equation  14b 

Typically  a  UCT  would  be  a  faint  object  that  appears  unexpectedly  in  the  sensor  field  of  view.  The  range  of 
Brightness  values  for  such  a  UCT  is  a  function  of  the  detection  limit  for  the  sensor.  This  range  is  typically  a  known 
value  and  is  randomly  distributed  between  the  range  of  Brightness  for  the  UCT  objects.  The  range  for  the 
resulting  values  for  changes  as  function  of  the  phase  angle.  These  values  for  are  used  to  compute  the 
empirical  PDF  for  /(  |  UCTj).  Note  that  this  PDF  also  becomes  a  function  of  the  phase  angle. 

With  the  above  calculations,  all  terms  on  the  right  hand  side  of  the  Bayes  expression  (Equation  8)  are  estimable  and 
considered  known.  This  allows  the  calculation  of  Brightness  evidence,  denoted  as  Xjk,  which  is  sent  to  the  GEO 
status  node.  Even  though  this  calculation  requires  the  determination  of  a  number  of  individual  terms,  it  is  a  minor 
computational  effort  for  a  desktop  computer.  It  is  represented  by  vertical  arrows  in  Eigure  6.1.  The  calculation 
corresponding  to  the  horizontal  arrows  is  described  in  Section  8.  In  this  regard,  note  that 

P{ANOMjk  I  rjk)  =  1  -  P{NOMjk  \  Equation  15 
=  1  -  Xjk 

8.  TRANSITION  PROBABILITY 

The  purpose  of  the  calculations  described  in  this  section  is  to  compute  the  terms  which  are  shown 

alongside  the  horizontal  arrows  in  the  Bayes  network  graph  (Eigure  6.1).  They  describe  the  probability  of  transition 
from  NOM  to  NOM  or  from  ANOM  to  NOM  in  the  interval  of  time  between  pass  k  and  pass  k+1  when  there  is  no 
data  collected.  This  interval  is  denoted  as  the  Inter-pass  Time  or  At/^in  the  following  description.  Accounting  for  the 
probability  of  change  of  state  during  the  Inter-pass  Time  is  necessary  for  the  following  reasons:  (1)  The  Inter-pass 
Time,  At/^  can  change  from  one  pass  to  the  next  by  over  a  magnitude.  Eor  example,  in  case  of  a  ground-based  sensor 
that  executes  an  east-west,  back-and-forth  sweeping  motion  to  perform  synoptic  search,  the  Brightness  Data  for  the 
satellites  located  at  high  zenith  angles  has  an  alternating  set  of  large  and  small  values  for  At/^;  (2)  There  is,  in  effect, 
no  collection  of  Brightness  Data  50%  of  the  time  due  to  the  daytime  gap;  (3)  The  At/^  is  routinely  significant  in 
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comparison  to  the  time  required  for  a  satellite  to  maneuver,  to  change  its  solar  panel  offset,  or  to  execute  a  1 80°  yaw 
turn  or  (i.e.  to  transition  towards  the  occurrence  of  an  ANOM  state  such  as  CrossTag,  OffsetChange  or  Turn).  The 
Atj^  can  also  be  significant  with  respect  to  a  duration  over  which  a  satellite  may  experience  degradation  of  its 
attitude  control,  progressing  towards  the  occurrence  of  the  ANOM  state.  Unstable;  (4)  The  Atj^  is  significant  in  order 
to  allow  different  UCT  objects  to  appear  in  the  successive  passes,  particularly  when  the  sensor  line-of-sight  can 
cross  a  debris  field  at  an  altitude  lower  than  the  geosynchronous  belt. 


For  the  purpose  of  this  section,  assume  that  Atj^  is  constant  between  all  night-time  passes  and  the  daytime  gap  is 
twelve  hours.  These  assumptions  are  unnecessary  in  the  actual  implementation  of  the  algorithm,  but  they  are  useful 
to  simplify  the  following  description.  Let  the  daytime  gap  be  denoted  as  AtQ^p- 


The  transition  probability  matrix  T  is  defined  using  the  following  matrix  equation  [References  15-17]. 


Where 


T  = 


Equation  16 


is  the  probability  of  persistence  of  NOM  state  (i.e.  no  change)  over  a  unit  interval  of  time 
is  the  probability  of  transition  of  ANOM  to  NOM  over  a  unit  interval  of  time 
is  the  probability  of  transition  of  NOM  to  ANOM  over  a  unit  interval  of  time 
is  the  probability  of  persistence  of  ANOM  over  a  unit  interval  of  time 


The  matrix  T  is  also  called  the  Markov  transition  matrix.  If  Atj^  =  1,  then 


N’T'k  N’T'k 

N^k+1  A^k+1 

A’T'k  A’T'k 

lN^k+1  A^k+l\ 


Equation  17a 


If  Atj^  ^  1,  the  Markov  transition  matrix  is  computed  using  the  Chapman-Kolmogorov  equation  [Reference  17]: 


N’T'k 

N^k+1 

A’T'k 

.N^k+1 


N’T'k 

A^k  +  1 

A’T'k 

A^k  +  1 


= 


Equation  17b 


The  (Atj^y^  power  of  the  Markov  matrix  T  is  computed  using  the  Cayley-Hamilton  theorem  [References  18-20]. 
This  comprises  a  three  step  procedure:  (1)  Compute  eigenvectors  and  eigenvalues  of  T;  (2)  Raise  the  eigenvalues  to 
the  (Atj^y^  power;  (3)  Compute  using  the  eigenvectors  for  matrix  T  and  its  eigenvalues  raised  to  the  (Atj^y^ 
power.  Note  that  each  row  of  the  Markov  matrix  T  sums  to  one  (see  Equations  6c  and  6d). 


A  salient  property  of  a  Markov  matrix  is  that  its  first  eigenvalue  is  equal  to  1  [References  18  and  21].  Also  the  sum 
of  two  eigenvalues  of  T  is  equal  to  the  sum  of  its  diagonal  elements  (i.e.  its  trace)  and  the  product  of  the  two 
eigenvalues  is  equal  to  the  determinant  of  T.  This  is  because  the  trace  and  the  determinant  of  a  matrix  are  its  first 
and  third  invariant,  respectively  [Reference  22].  Thus,  the  two  eigenvalues  for  matrix  T,  namely  A^and  A2  are  given 
directly  by  Equation  1 8  and  two  eigenvectors  are  computed  by  solving  two  sets  of  2x2  simultaneous  equations : 


=  1  Equation  18a 

^2  =  Equation  18b 

Using  the  same  logic  to  the  transition  probability  for  the  daytime  gap  by  replacing  the  (Atj^y^  power  of  the  Markov 
matrix  T  to  the  power: 


N’T  N’T 

N^Gap  A^Gap 

A’t  A’t 

N^Gap  A^ Gap 


Equation  19 


The  matrix  elements  on  the  left  hand  side  of  Equation  19  are  defined  in  analogous  manner  as  Equation  17b;  except 
that  the  right  superscript  and  subscript  is  replaced  by  the  descriptor  Gap  instead  of  the  pass  numbers  k  and  Note 
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that  is  user  input,  as  defined  in  Section  4.  This  allows  the  calculation  of  A^cap  using  Equation  6c.  However, 

the  matrix  elements  in  the  second  row  of  Equation  19  are  unknown.  They  are  solved  using  the  user  input  defined  in 
Section  4  for  the  long  run  probability  that  the  satellite  is  NOM  at  any  time.  This  is  described  next. 


In  the  limit,  long  term  is  defined  as  being  at  least  two  orders  of  magnitude  longer  thanAt^^p,  which  is  3  to  6 
months.  The  long  run  probabilities  may  now  be  visualized  as  saturated  or  steady-state  values  resulting  from 
successive  multiplication  of  daytime  gap  transition  probability  matrix  ad  infinitum: 


Un  (1  — 

.(1  —  TT^)  TTyi 


N'T'  N'T' 

N^Gap  A^Gap 

A’p  A’p 

N^Gap  A^ Gap 


Equation  20 


Note  that  is  the  user  input  value  of  long  run  probability.  Note  that  the  rows  of  both  matrices  are  consistent  with 
Equation  6.  However,  in  the  steady  state,  is  also  equal  to  (1  -  tt^)  and  thus  Equation  20  is  singular.  It  is  however, 
equivalent  to  a  left  eigenvector  form  as  follows: 


[Un  =[n:M 


t^a] 


'N’T' 

N^Gap 

A’p 

N^Gap 


N'p 

A^  Gap 

A’p 

A^  Gap 


Equation  21 


Equation  21  is  solved  using  the  user-input  values  for  and  N^^cap  ^  determine  every  term  in  the  equation.  Thus, 
this  equation  is  the  starting  point  for  the  calculation  of  the  transition  probability  matrix  T  as  the  (l/Atoap)*  root  of  the 
transition  probability  matrix  for  the  daytime  gap.  This  is  performed  using  the  Cayley-Hamilton  theorem. 


Note  how  the  procedure  is  generalized  when  Atj^  and  AtQ^p  are  non-uniform.  Even  though  this  calculation  requires 
the  determination  of  a  number  of  individual  terms,  it  is  a  minor  computational  effort  for  a  desktop  computer  because 
it  involves  a  set  closed-form  analytical  procedure;  namely  the  use  of  Chapman-Kolmogorov  equation  to  calculate 
the  multi-step  transition  matrices  and  the  use  of  Cayley-Hamilton  theorem  to  perform  this  calculation  quickly. 

9.  MARKOV  CHAIN 


Markov  chain  is  utilized  to  proceed  from  one  GEO  status  node  to  the  next  by  combining  the  calculations  from  the 
application  of  Bayes  theorem  and  the  transition  probability.  The  two  calculations  taken  together  allow  computation 
of  the  term  )  using  the  assumption  that  the  NOM  or  ANOM  state  for  a  satellite  during  pass  k-\-l 

depends  only  on  its  state  during  pass  k  (Section  4).  This  is  Markov  chain  of  memory  1,  or  order  1,  or  simply  a  first 
order  Markov  chain  [References  15  and  17].  It  is  possible  to  utilize  a  higher  order  Markov  chain,  but  it  needs  joint 
probability  distribution  functions.  This  requires  larger  amounts  of  data  as  input,  which  conflicts  with  the  user  need  to 
produce  credible  statistical  assessment  while  minimizing  the  input  data  required.  The  first  order  Markov  chain 
provides  matrix  equations  for  the  calculation  of  P(  )  and  for  the  calculation  of  the  numerator  in  Equation 

9  for  pass  (k-i-1)  as  follows: 


[P(  )  P(  )]  =  [Xjk 


N’T’k 

N^k+1 

A’pk 

.N^k+1 


N’T’k 

A^k+1 

A'j'k 

A^k+1 


Equation  22a 


/(  I  )  P{  = 

/(  I  XXj,  +  (  1  -  Xj^)  )  Equation  22b 

All  entities  on  the  right  hand  side  of  Equation  22  are  known.  This  allows  the  calculation  of  P(  ),  which  is 

used  in  the  application  of  Bayes  theorem  as  per  Equation  8.  The  Bayes  network  is  begun  with  a  known  state  for  a 
satellite  at  /c  <  0  and  the  Markov  chain  calculation  allow  the  calculations  to  traverse  the  network  forward.  Its  result 
is  to  produce  a  vector,  denoted  as  Xy,  which  contains  a  sequence  of  numerical  values  of  P  (  NOMjj^  \  )  or  Xjk  for 

all  k,  both  in  the  prior  data  and  the  new  data  (Eigure  3.1).  Or,  the  length  of  vector  Xj  is  equal  to  the  running  total 
number  of  Brightness  Data  points  in  the  prior  data  and  new  data.  The  belief  or  marginal  probability  of  the  NOM 
state  is  the  expected  value  for  P  (  NOMjj^  I  ^7c  )•  Or, 
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Belief  (NOMjk)  =  IEi?[  P  (  NOMji,  |  /?  =  [ Xj  ] 


Equation  23 


Where  denotes  the  expected  value  (i.e.  mean)  with  respect  to  all  values  of  R.  The  belief  can  be  computed  for  the 
prior  data  and  compared  with  the  new  data  (see  Section  15).  Or  it  can  be  computed  for  a  moving  set  of  passes  (e.g. 
for  a  moving  block  of  thirty  passes  such  as  from  0  <  /c  <  29,  1  <  /c  <  30,  etc. 


10.  OUTLIERS  AND  ANOMALY  DETECTION 

If  the  belief  changes  from  the  prior  data  to  the  new  data,  it  implies  anomaly  or  change  in  the  satellite  (see  Equation 
23).  However,  the  data  set  corresponding  to  k  >  0  can  be  very  small.  Eor  example,  in  case  prior  data  that  spans  two 
weeks,  there  would  be  about  30  points  of  Brightness  Data  for  each  of  the  three  regions  in  Eigure  2.2,  or  the  length  of 
vector  Xj  for  each  region  of  the  prior  data  would  also  be  30.  The  new  data  for  Xjk  for  each  of  the  three  regions  may 
be  collected  at  a  rate  of  three  points  per  day.  Thus,  challenge  is  to  identify  change  in  a  satellite  based  on  reference 
data  set  that  has  about  30  points  and  new  data  that  is  as  small  as  three  points.  In  such  cases,  the  use  of  Equation  23  to 
calculate  belief  is  inadequate.  As  a  result,  the  anomaly  detection  is  performed  using  a  combination  of  three  types  of 
evidence;  namely  Xj^  and  two  likelihood  ratios  determined  using  cluster-based  evidence  and  model-based  evidence. 

The  first  evidence  is  to  determine  if  the  Xjk  for  the  current  pass  is  an  outlier  with  respect  to  the  prior  data  for  Xjk,  k  < 
0.  This  evaluation  is  performed  separately  for  the  specular  and  each  Lambertian  region  so  that  the  outliers 
corresponding  to  the  various  anomalies  described  in  Section  5  can  be  identified.  This  is  described  in  Section  11. 
The  second  evidence  is  based  on  the  values  of  order  to  identify  the  ripple  effect  of  the  outlier  in  Xj^.  This  is 

described  in  Section  12.  The  third  evidence  is  to  perform  model  based  correction  and  to  re-evaluate  Xj^  to  assess  if  it 
continues  to  be  an  outlier.  Eor  example,  in  case  of  CrossTag,  the  term  in  Equation  2b  can  be  temporarily 
replaced  by  the  predicted  brightness  for  the  cluster  peer  satellite.  This  is  described  in  Section  13. 

These  three  sets  of  evidence  are  combined  and  propagated  as  per  the  Bayes  network.  This  is  described  in  Section  14. 
A  change  in  satellite  behavior  is  detected  if  the  evidence  of  an  outlier  persists  over  multiple  cycles.  In  other  words,  it 
is  an  iterative  assessment.  This  calculation  can  be  performed  even  if  the  new  data  comprises  very  few  passes.  These 
calculations  are  simpler  but  they  do  not  allow  well -quantified  hypothesis  tests,  which  are  essential  in  order  to 
generate  probabilities  of  Type  I  and  Type  II  errors,  a  and  p,  which  are  related  to  the  confidence  and  power  of  the 
test.  Eor  example,  a  is  the  value  of  the  probability  that  a  satellite  would  be  deemed  ANOM  when  it  actually  NOM. 
This  is  described  further  in  Section  16. 

II.  USE  OF  Z-SCORE  FOR  DETECTION  OF  OUTLIERS 

The  rationale  for  this  test  is  to  consider  the  (l-Xjk)  values  for  the  prior  data  to  be  the  status  quo  and  examine  the  (1- 
Xjk)  for  each  pass  k  >  0  in  the  new  data  to  assess  if  there  has  been  a  change  in  (l-Xjk)  from  the  status  quo.  The 
status  quo  is  defined  through  the  mean  and  standard  deviation  of  (1-  Xjk)  for  passes  k  <  0  for  the  prior  data. 
Specifically: 


^k<o  ^y/c)] 


k  <  0  Equation  24 


(4<of  =  [(1  -  Xj,)] 


k  <  0  Equation  25 


AMOS  2014  Technical  Conference 


Where 


Ej^  denotes  the  expected  value 

is  the  mean  of  (l-Xjk)  for  the  prior  data,  k  <  0 

is  the  standard  deviation  of  (l-Xjk)  for  the  prior  data,  k  <  0 

The  first  step  is  to  estimate  the  quantities  and  based  on  a  limited  set  of  Xj^  values  for  the  prior  data.  At 
present,  this  calculation  is  performed  by  assuming  that  the  mean  and  variance  of  Xj^  are  constant  for  all  k  <  0  but 
depend  on  j .  These  estimates  for  the  mean  and  variance  are  denoted  as  and  ,  respectively. 

In  order  to  determine  whether  or  not  there  has  been  a  change  in  Xjk  from  the  status  quo,  a  basic  z-score  measure 
is  defined  in  the  following  [Reference  23].  It  is  useful  when  only  a  few  observations  of  new  data  are  available.  The 
magnitude  of  this  z-score  is  meant  to  prompt  the  user  to  assess  if  the  satellite  has  changed  relative  to  the  status  quo. 
Specifically,  the  z-score  for  pass  k  is  given  by: 


Cfc  =  j  Equation  26 

^0 

The  user  chooses  a  threshold  level  for  the  magnitude  of  outside  of  which  the  satellite  is  expected  to  be  ANOM. 
This  calculation  is  performed  at  each  pass  k  and  the  result  is  a  time  evolution  of  the  z-score  based  on  which  a  user 
may  make  a  judgment  in  regards  to  a  change  in  a  satellite.  Note  that  for  this  computation,  the  user  can  choose  to 
update  the  values  of  mean  and  variance  that  are  used  as  inputs  to  the  z-score  calculation. 

Note  that  the  knowledge  of  z-score  evolution  over  time  is  useful  for  a  basic  assessment  of  the  anomaly  type  for  the 
generation  of  a  message  for  use  in  the  Bayes  belief  propagation.  For  example,  if  the  z-score  outliers  are  only  for  the 
passes  corresponding  to  observation  geometry  for  the  specular  region,  it  is  likely  to  be  caused  by  OffsetChange.  If 
the  outliers  occur  sporadically  within  a  night,  it  may  be  due  to  repeated  occurrences  of  CrossT ag .  If  an  outlier  only 
occurs  for  an  isolated  data  point  and  then  corrects  itself,  it  may  be  due  to  a  mistag  caused  hy  UCT.  If  the  outliers  are 
persistent  at  all  phase  angles  and  possess  a  structure,  it  may  be  due  to  a  180°  yaw  Turn  of  the  satellilte.  If  the 
outliers  are  persistent  at  all  phase  angles  and  lack  a  structure,  it  may  be  because  the  satellite  is  Unstable,  etc. 

12.  CLUSTER  BASED  SUPPORTING  EVIDENCE  OF  ANOMALY 


As  illustrated  in  Figure  2.1,  the  Brightness  Data  can  be  collected  for  all  satellites  in  a  cluster  simultaneously. 
Accordingly,  it  is  assumed  that  the  Brightness  Data  for  the  entire  cluster  is  available  to  the  statistical  assessment 
algorithm  for  each  pass  (Figure  6.1).  Indeed,  the  statistical  assessment  will  run  in  a  sequence  for  one  satellite  in  the 
cluster  after  another,  creating  a  separate  chain  of  data  for  each  satellite  (Figure  12.1a).  The  calculation  of  vector  Xj 
and  the  z-score  values  can  thus  be  performed  sequentially  for  each  satellite  Cj  that  belongs  to  the  cluster.  This  is 
practical  because  of  the  sparse  rate  at  which  the  Brightness  Data  is  collected.  The  duration  of  time  between  passes  is 
expected  to  be  lOx  longer  than  the  time  required  to  compute  Xjk  and  values  for  each  satellite  in  the  cluster. 


The  cluster  data  for  Xjk  and  z-score  provides  a  second  opportunity  to  gather  evidence.  For  example,  if  the  z-score  is 
an  outlier  for  all  satellites  in  the  cluster,  it  may  be  due  to  a  noisy  observation  or  due  to  the  crossing  of  a  debris  field 
(UCT)  in  the  line  of  sight  of  the  sensor.  If  the  outliers  occur  simultaneously  for  cluster  peers,  it  may  be  due  to 
CrossTag.  Such  cluster-based  evidence  is  expressed  in  terms  of  a  likelihood  ratio  such  as: 


Cluster  ^  SdCrossTag)  ^  p{aN0Mj„\  r \ 
£  (no  CrossTag)  P(W0Mjfc|  | 

P(ANOMj,\  rj„)  =  1  -  P{N0Mj„\  rj,)  =  1  -  Xj, 


Equation  27  a 


Equation  27b 


Where  £  (^CrossTag)  is  likelihood  for  the  presence  of  CrossTag,  and  is  the  message  sent  to  the  GEO 

status  node  for  satellite  Cj  pass  k.  Use  of  during  the  belief  propagation  is  described  in  Section  14.  Note  that 

multiple  such  messages  can  be  derived  from  the  cluster-based  evidence  and  sent  to  the  GEO  status  node. 
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Figure  12.1a:  Bayes  network  for  each  satellite  in  the  cluster  are  independent  of  each  other 


Figure  12.1b:  Predictive  model  can  be  applied  to  each  satellite  in  the  cluster 

13.  MODEL  BASED  SUPPORTING  EVIDENCE  OF  ANOMALY 

The  present  work  utilizes  the  two -facet  model  to  perform  inversion  of  prior  data  to  determine  the  pose-dependnent 
body  and  solar  panel  albedo-area  products  for  three-axis  stabilized  GEO  satellites  (Figure  12.1b).  The  Brightness 
prediction  for  the  calculation  of  rjk  is  performed  by  utilizing  the  previously  computed  albedo -area  products  along 
with  the  specific  observation  conditions  for  pass  k  in  the  new  data.  This  calculation  is  a  generic  process  and  it  is 
applied  sequentially  to  each  satellite  in  the  cluster.  The  ability  to  use  the  two-facet  model  on  demand  provides  a 
third  opportunity  to  gather  evidence,  which  can  be  communicated  back  to  the  GEO  status  node  as  a  message.  This 
evidence  is  extracted  as  a  likelihood  ratio  based  on  the  PDF  of  the  Brightness  Ratio  as  follows: 
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For  example,  consider  a  situation  when  the  cluster  based  evidence  indicates  possibility  of  CrossT ag  anomaly  for 
satellites  Cj  and  •  If  the  CrossT  ag  was  true,  the  Brightness  Ratio  for  the  two  satellites  utilized  in  the 
calculation  of  Xij  and  have  been  calculated  using  incorrect  values  of  Vjj^  and  as  follows  (Equation  2): 

For  satellite  Cj :  Vjj^  =  and  rjj^  =  1  —  Vjj^ 

^mjk 

For  satellite  and  =  1  — 

The  Predictive  Model  may  be  used  to  perform  a  numerical  trial.  If  there  was  CrossT  ag,  a  trial  for  the  resolution  of  a 
crosstag  can  be  performed  by  correcting  the  values  of  Vjj^  and  .  This  requires  that  predicted  values  of 

Brightness  used  in  denominator  for  the  calculation  of  Vjj^  and  are  switched  from  satellite  Cj  to  and  vice 

versa.  Then, 


For  satellite  Cr.  Vju  =  For  satellite  : 

^  Imjk  J 


jPeer 

^jPeer  _  ojk 
^jk  jPeer 

^mjk 


The  likelihood  ratio  for  CrossT  ag  may  now  be  computed  using  the  PDF  for  the  corrected  values  of  the  Brightness 
Ratio  for  the  cluster  peers  as  follows: 


Model  ^  CrossTag  \  corrected  rjj^  and  rj^^^  ) 

£  (CrossTag  \  corrected  rjj^  and  rj^^^  ) 

P  {corrected  rjj^  \  NOMjj^)  P(corrected  rj^^^  \ 
p{corrected  rjj^  \  ANOMjj^)  P(corrected  \  ANOM^^^'^) 


Equation  28 


Where  is  the  message  from  the  model  to  the  GEO  status  node  for  satellite  Q  pass  k.  Use  of  during 

the  belief  propagation  is  described  in  Section  14.  Note  that  multiple  such  messages  can  be  derived  from  the  model - 
based  evidence  and  sent  to  the  GEO  status  node. 


14.  BELIEF  PROPAGATION 

Belief  is  marginal  probability.  Accordingly,  Equation  23  can  be  used  to  determine  belief  at  each  GEO  status  node. 
However  a  direct  use  of  Equation  23  requires  a  minimum  of  30  points  of  new  data.  In  order  to  be  able  to  assess 
belief  even  when  the  size  of  new  data  is  small,  the  procedure  described  in  this  section  is  utilized.  It  propagates  a 
representative  measure  of  belief  by  combining  the  three  different  values  of  evidence  described  in  Sections  11-13. 
Specifically:  (1)  z-score  for  P(NOMjj^  |  rj^')  ,  which  identifies  outliers.  These  are  the  passes  where 
P{ANOMjj^  I  rjj^)  is  high  ;  (2)  Message  ,  which  is  a  cluster-based  likelihood  ratio  for  ANOMjj^;  (3)  Message 

rrip^^dei^  which  is  a  model-based  likelihood  ratio  for  the  resolution  of  ANOMjj^.  Figure  14.1  illustrates  a  process  for 
the  propagation  of  representative  value  for  belief  by  combining  the  three  different  values  of  evidence  values  at  each 
GEO  status  node  as  follows  [Reference  5]: 

(ANOMj^)  =  Equation  29 

Where  is  a  representative  measure  of  belief  or  score  for  ANOM  at  pass  k.  A  user  can  define  threshold  values  for 
iB  .  The  values  of  iB  cross  can  be  monitored  in  real-time  and  a  warning  issued  if  the  threshold  is  crossed  repeatedly. 
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^Model 

rrik 


^Cluster 


15.  DECISION  TO  MOVE  THE  TIME  SLIDER 


At  the  root  of  the  mathematics  presented  in  the  preceding  sections  is  the  prior  data,  which  is  used  to  calculate 
/( ^jk  I  NOMj).  Since  the  prior  data  is  for  the  NOM  state,  which  itself  can  evolve  over  time,  it  is  necessary  to 
exercise  due  diligence  prior  to  denoting  that  the  satellite  is  NOM  for  a  set  of  passes  (i.e.,  over  an  interval  of  time).  It 
is  further  necessary  that  this  set  comprise  at  least  30  passes  as  described  in  the  following. 

It  is  important  to  attain  a  balance  between  the  probabilities  of  false  alarm  and  the  probability  of  missing  anomalous 
behavior  in  a  way  that  fulfills  the  user  requirements.  These  two  values  of  probability  cannot  be  obtained  using  the  z- 
score  approach  described  in  Section  11.  These  two  probabilities  can  be  determined  by  developing  suitable 
hypothesis  tests  and  computing  the  estimates  of  Type  I  and  Type  II  errors  [Reference  23].  In  order  to  compute  these 
probabilities,  it  is  necessary  to  develop  certain  distributional  assumptions  in  regards  to  for  the  prior  data  (k  <  0), 
and  employ  a  z-score  measure  for  the  new  data  (k>0)  for  which  a  probability  distribution  can  be  obtained  or 
approximated.  This  requires  sufficiently  large  number  of  points  for  Xj^  for  the  prior  data  and  new  data  (such  as  a 
sample  size  of  30  or  more  for  each).  Upon  collecting  such  data  and  determining  that  the  satellite  behavior  is  NOM, 
the  time  slider  can  be  moved  forward  to  new  position  up  to  which  the  satellite  is  shown  to  be  NOM. 

The  calculation  of  probability  distribution  for  the  z-score  would  be  easier  if  Xj^  were  independent  variables. 
However,  this  is  not  likely  to  be  the  case  as  each  Xj^  is  computed  using  input  from  Thus,  the  present 

problem  is,  in  fact,  a  part  of  the  eternal  statistical  conflict  between  small  sample  sizes  and  margin  of  error. 

Section  11  defines  the  z-score  measure  using  Equation  26.  Its  simplicity  comes  at  the  expense  of  not  being  able  to 
quantify  the  probabilities  of  Type  I  and  Type  II  errors.  To  this  end,  a  second  z-score  measure  77^  is  defined  as 
follows  [Reference  26] : 

/  xJ  - 

Vk  =  J — Equation  30 

^k<0  /  _ 

/  V/c+1 

where  denotes  the  sample  mean  of  observations  Xjo,...,Xjk,,  k  >0.  If  the  Xjo,...,Xjk  were  independent  and 
identically  distributed,  this  statistic  would  have  an  approximate  standard  normal  distribution  for  large  enough  k  and 
could  be  used  as  a  test  statistic  for  a  hypothesis  test.  This  test  would  allow  the  calculation  of  Type  I  and  Type  II 
errors  using  standard  theory.  Specifically: 

Ho  :  Mfc>o  is  equal  to  Ha :  mean  Mfc>o  is  not  equal  to 

where  Mfc>o  “  ^k[^jk>]  ,k>0  Equation  31 

However,  the  values  for  Xjo,...,Xjk  are  not  independent.  They  are  more  likely  to  be  a-mixing  [Reference  26].  The 
assumption  of  a-mixing  indicates  an  approximate  independence.  It  is  reasonable  to  assume  that  the  sequence  Xjk  is 
m-dependent,  which  is  defined  as  follows: 
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(Xjo,. .  .,Xjh  )  is  independent  of  (Xjn,  Xj(n+i), . )  whenever  n-h  >  m,  for  some  integer  m  and  any  h. 

If  so,  the  a-mixing  property  follows  accordingly.  Furthermore,  assuming  that  Xjk  is  stationary  (i.e.,  the  joint 
probability  distribution  of  (Xjn,.---,Xj(n+k)  )  does  not  depend  on  n),  a  Central  Limit-type  theorem  for  stationary  a- 
mixing  sequences  can  be  employed  so  that  77^  has  an  approximately  normal  distribution,  with  replaced  by  as 
follows  [Reference  26] : 

limitj^^oo  +  ^71+^72  +  — f  ^jk)  —  (  Equation  32a 

where  (  =  IE  [(Xyo  - 

+  2  lE[(A’yo  -  Equation  32b 

Equation  32  is  an  adjustment  of  the  variance  of  each  Xji  to  account  for  the  dependence.  It  is  necessary  to  account  for 
this  dependence  because  it  affects  the  accuracy  of  Type  I  and  Type  II  error  estimates,  compromising  the  calculation 
of  the  probability  for  the  false  positives  and  false  negatives. 

Since  77^  is  approximately  normally  distributed  for  large  enough  k,  the  hypothesis  is  rejected  when  the  observed  test 
statistic  value  is  such  that 

vi  <  -Za/2  or  r]l  >  z„/2 

where  Zc^/2  denotes  the  upper  (100)(a/2)*  normal  percentile.  Rejecting  the  null  hypothesis  means  that  the  test 
procedure  concludes  that  the  mean  of  the  Xjk  for  passes  k  >  0  is  different  from  its  value  for  passes  k  <  0.  When  the 
hypothesis  test  is  performed  at  a  significance  level  a,  there  is  a  probability  of  a  that  the  test  procedure  will  cause  the 
incorrect  assertion  that  the  mean  Xjk  for  new  data  (  k  >  0  )  is  not  equal  to  its  estimated  mean  based  on  the  prior 
data,  Mfc<o-  This  is  the  Type  I  Error. 

To  quantify  a  Type  II  Error,  suppose  that  the  user  requirement  is  to  attain  at  least  (1-p)  probability  that  the  chosen 
test  procedure  correctly  rejects  the  null  hypothesis,  when  the  true  mean  based  on  the  process  for  passes  k  >  0  is  not 
q{^<o  but  rather  it  is  +  6.  In  this  case,  it  is  necessary  to  consider  the  sample  mean  of  (Xjo+. .  .+Xjk)/(k-Fl),  where 
sample  size  k-^I  is  at  least  (za/2+zp)^ajV(  5^  )  [Reference  26].  The  justification  of  this  power,  as  well  as  the 
distribution  of  the  test  statistic  itself,  assumes  that  the  estimates  of  and  of  are  "known"  quantities. 

There  are  two  drawbacks  of  this  method.  First,  the  test  is  predicated  on  full  confidence  in  the  estimates  and 
The  margins  of  error  associated  with  these  estimates  play  an  underlying  role  in  the  error  rates  of  the  test. 
Second,  estimating  is  a  complex  problem  related  to  estimating  a  covariance  matrix  of  a  multivariate 

distribution  and  it  is  difficult  to  compute  when  the  data  set  is  small.  As  a  practical  matter,  if  the  prior  data  set 
contains  too  few  observations,  this  complicates  the  estimation  of  and/or  increases  the  margin  of  error  for  the 
estimate.  It  is  anticipated  that  this  situation  will  lend  itself  to  the  bootstrap  resampling  technique  in  order  to  evaluate 
and  quantify  confidence  in  this  estimate  [Reference  24] . 
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16.  NEAR-REAL  TIME  ASSESSMENT 


Figure  16.1  illustrates  how  the  Control  Logic  in  the  Statistics  Model  is  posed  such  that  it  is  feasible  to  provide  a 
near-real  time  assessment  of  NOM  or  ANOM  states  for  GEO  satellites  as  each  new  point  of  Brightness  Data  is 
received.  The  procedure  is  near-real  time  in  that  the  assessment  can  be  obtained  well  in  advance  of  the  time  when 
the  next  data  point  for  Brightness  Data  is  expected  to  become  available.  The  salient  components  of  the  Statistics 
Model  are  identified  with  box  numbers  in  the  flow  chart.  These  boxes  can  be  mapped  to  various  figures  and 
equations  in  the  paper,  as  summarized  in  Table  16.1.  Almost  all  calculations  comprise  small  computational  effort  for 
a  desktop  computer.  The  most  computationally  intensive  part  of  the  method  is  the  utilization  of  the  two -facet  model 
for  inversion  of  the  prior  data  (Figure  3.1),  taking  several  minutes  per  satellite.  However,  this  calculation  is  required 
only  when  the  time  slider  is  moved  to  a  new  location. 

The  utilization  of  a  Bayes  network  for  statistical  assessment  has  important  collateral  benefits  as  follows.  It  generates 
a  chronological  chain  for  photometric  characterization  data  for  a  satellite,  which  lends  itself  naturally  to  a  “last-in 
first-out”  type  of  data  archival  structure.  A  chain  of  such  data  can  be  generated  for  each  satellite  of  interest  and  then 
used  for  its  life  cycle  assessment.  Furthermore,  the  Bayes  network  structure  for  photometry -based  assessment  can  be 
combined  with  the  Bayes  network  structure  for  metrics-based  assessment  in  order  to  create  the  ability  to  perform  a 
metric-photometric  assessment  of  satellite  status  in  near-real  time  [Reference  4]. 
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Figure  16.1:  Control  Logic  in  the  Statistics  Model 
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Table  16.1:  Mapping  of  the  flow-chart  boxes  to  technical  description 


Box 

Title 

Mapping  to  the  technical  description 

1 

Initialize  to 

Section  3,  Timer  slide  definition 

2 

Set  pass  k  =  0 

Figure  3.1 :  Prior  data  for  k  <  0.  New  data  for  k  >0 

3 

Compute  prior  distribution 

Sections  5  and  7:  PDF  for  NOM  and  ANOM 

4 

Compute  transition  probability 

Section  8:  Chapman-Kolmogorov  equation 

5 

Receive  data  for  new  pass 

Section  2:  Character  of  Brightness  Data 

6 

Compute  P(NOMk  1  rk) 

Sections  6,7  and  9:  Bayes  network,  Bayes  theorem  and  Markov  Chain 

7 

Determine  belief 

Sections  10-13:  Outlier  detection,  z-score,  cluster  and  model  evidence 

8 

Notify  user 

Section  14:  Propagate  belief  and  notify  user 

8 

Update  to 

Section  15,  Movement  of  the  time  slider 

17.  ONGOING  WORK  AND  CLOSURE 

In  light  of  the  formulation  for  belief  propagation  and  anomaly  detection  described  in  the  preceding  sections,  future 
work  entails  computing  the  probabilities  of  Type  I  and  Type  II  errors  in  light  of  the  unavoidable  small  sample  sizes. 
Even  if  a  large  enough  training  set  of  prior  data  exists  for  belief  in  satellite  status,  the  user  does  not  want  to  wait  for 
a  large  number  of  passes  in  order  to  assess  whether  the  satellite  has  changed  from  a  NOM  state. 

To  this  end,  the  future  work  will  attempt  to  address  this  problem  in  the  context  by  exploring  how  the  margin  of  error 
can  be  mitigated  in  two  ways.  First,  bootstrap  analysis  will  be  used  to  improve  the  estimate  of  distributional 
properties  of  test  statistics  (77^).  Once  the  distributional  properties  are  established,  the  calculation  for  the  margin  of 
error  can  be  improved  in  order  to  provide  the  user  with  the  critical  values  associated  with  rejecting  steady- state. 

The  use  of  an  alternate  functional  (i.e.,  a  function  whose  domain  is  a  set  of  functions)  with  more  desirable  statistical 
properties  will  be  explored  in  order  to  represent  the  belief  process.  Here  “more  desirable”  means  that  the 
distributional  properties  of  a  test  statistic  derived  from  this  functional  can  be  calculated  within  the  constraints  of  the 
current  problem  of  anomaly  detection.  For  example,  such  a  statistic  may  exploit  potential  independence  of 
increments  (such  as  Xj^  —  ^;(/c+i)),  independence  of  time  blocks  (such  as  belief  gathered  at  non -overlapping  time 
frames),  or  other  measures  of  stability  (e.g.,  sample  median). 

The  distribution  of  the  belief  process  in  the  larger  context  of  stochastic  processes  will  be  explored.  Statistical 
analysis  involving  model  fitting  is  the  first  step  toward  this  goal,  including  autoregressive  or  other  mean  reverting 
models.  Software  already  exists  to  obtain  parameter  estimates  for  such  processes,  and  as  such  a  comparison  of 
estimated  parameters  in  the  prior  with  those  of  new  data  could  lend  itself  to  improvement  in  the  definition  of  the 
PDF  for  prior  and  new  data  and  an  additional  test  for  change  in  belief  in  the  status  of  a  satellite. 
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APPENDIX  A:  NOTATION 


Text  notation: 


Brightness:  Value  of  Brightness  Data  at  a  single  data  point 

Brightness  Data:  Single  point  brightness  data  point  collected  by  a  ground  or  space-based  sensor  during  the  routine 
synoptic  search  operation.  This  data  is  collected  along  with  the  angles  only  metric  data 
Brightness  Ratio:  Ratio  of  observed  Brightness  to  predicted  Brightness 
Body:  Body  of  a  GEO  satellite 

Change:  Modification  of  state  of  a  satellite  from  NOM  to  ANOM  or  vice  versa. 

Change  Detection:  To  recognize  that  the  state  of  satellite  has  undergone  Change 
Cluster  Peer:  A  pair  of  satellite  in  a  cluster  can  be  cross-tagged 
Inter-pass  Time:  Temporal  spacing  between  successive  orbital  passes 

NOM:  Nominal  status  of  a  satellite.  This  is  when  the  correlation  coefficient  between  the  observed  Brightness  and 
the  predicted  Brightness  exceeds  a  user-defined  threshold  limit. 

ANOM:  Anomalous  status  of  a  satellite.  This  is  when  the  correlation  coefficient  between  the  observed  Brightness 
and  the  predicted  Brightness  is  below  the  user-defined  threshold  limit  for  NOM. 

Panel:  Panel  term  for  a  GEO  satellite  (it  combines  the  effect  of  both  solar  panels  into  a  single  term) 

PDE:  Probability  distribution  function 

Signature  Data:  A  sequence  of  Brightness  measurements  collected  by  a  dedicated  sensor  during  a  single  pass  for  a 
target  satellite.  Eor  GEO  satellites,  such  data  is  collected  at  a  frame  rate  such  as  one  data  point  per  minute,  etc. 

Mathematics  notation: 

®  is  a  representative  measure  of  belief  or  a  score  for  ANOM  at  pass  k 
Cj  =  jth  satellite  in  a  GEO  cluster 

£  (CrossTag)  is  likelihood  for  the  presence  of  CrossTag 
Ej^  denotes  the  expected  value 

denotes  the  expected  value  (i.e.  mean)  with  respect  to  all  values  of  R 
Ey  =  Expected  value  or  mean  for  a  quantity  irrespective  of  y 
/( 7}/^)  is  the  PDE  for  the  Brightness  Ratio  for  pass  k  for  satellite  Cj 

/( rjk  I  ANOMj)  is  the  PDE  for  the  Brightness  for  pass  k  Ratio  given  the  satellite  Cj  is  ANOM 
/( rjk  I  CrossTag j)  is  the  PDE  of  Brightness  Ratio  when  Cj  is  cross-tagged  with  a  peer 
/( rjk  I  NOMj  )  is  the  PDE  for  the  Brightness  Ratio  for  pass  k  given  the  satellite  Cj  is  NOM, 

/( rjk  I  OffsetChangej)  is  the  PDE  of  Brightness  Ratio  when  the  Cj  panel  offset  is  modified 
/( Ik  I  Turrij)  is  the  PDE  of  Brightness  Ratio  when  Cj  is  turned  180° 

/( rjk  I  UCTj)  is  the  PDE  of  Brightness  Ratio  when  mistag  is  caused  by  an  unexpected  UCT 
/( ^jk  I  Unstablej)  is  the  PDE  of  Brightness  Ratio  when  Cj  is  unstable 
Imjk  is  the  Brightness  predicted  by  the  model  for  satellite  j  and  pass  k 
lojj^  is  observed  brightness  for  satellite  j  and  pass  k 
j  is  the  index  for  satellite  number  in  a  GEO  cluster 

or  So  is  the  standard  deviation  of  Xjk  for  the  prior  data,  k  <  0 
is  standard  z-score  for  anomaly  detection 
gl  is  another  z-score  measure  utilized  to  forward  the  time  slider 

k  =  Index  for  an  orbital  pass  number.  The  time  slider  origin  is  k  =  0.  Prior  data  is  for  k  <  0.  New  data  is  for  k  >  0. 
^Cluster  j^essage  based  on  cluster  based  evidence 
^M^odei  jYiessage  based  on  model  based  evidence 
£  {CrossTag)is  likelihood  for  the  presence  of  CrossTag 
The  long  run  probability  that  a  satellite  is  NOM  at  any  time 
P{Crosstagj)  is  the  long  run  probability  of  cross-tag  for  Cj 

P(ANOMj)  is  the  long  run  probability  that  the  satellite  is  ANOM,  also  denoted  as  P(  ANOMj')  =  1  - 
P(  NOMj)  is  the  long  run  probability  that  the  satellite  is  NOM,  denoted  as  tt^,  and 
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P  (  NOM^j  I  )  =  Probability  that  satellite  Ci  is  NOM  after  pass  k  given  the  Brightness  Ratio  rkj 
P{OffsetChangej)  is  the  long  run  probability  of  change  in  the  panel  offset  angle  for  satellite  Q 
P(  Turrij)  is  the  long  run  probability  of  a  180°  turn  for  Q 

P(  UCTj)  is  the  long  run  probability  of  mistag  for  Q  due  to  the  appearance  of  UCT  in  the  sensor  FOV 
P(  Unstablej)  is  the  long  run  probability  that  satellite  Q  is  unstable 
P(X  =  x)  =  Probability  that  an  outcome  for  a  random  variable  X  equals  v 
P(X  =  X  \  Y  =  y)  =  Conditional  probability  for  X  =  v  given  Y=y 

rjk  =  Ratio  of  observed  Brightness  of  satellite  j  and  pass  k  to  predicted  Brightness  of  satellite  j  and  pass  k.  This  ratio 
is  defined  only  at  the  orbital  location  when  Brightness  data  is  collected 
T  is  the  pass-to-pass  transition  probability  matrix 

A^k+i  —  ^  ~  N^k+i^  aT  is  the  probability  of  persistence  of  ANOM  over  a  unit  interval  of  time 
N^k+i  is  the  probability  of  transition  of  ANOM  to  NOM  over  a  unit  interval  of  time 
A^k+i  is  the  probability  of  transition  of  NOM  to  ANOM  over  a  unit  interval  of  time 
Mcap  is  the  daytime  gap 

Atk  is  the  temporal  spacing  between  pass  k  and  pass  (k+1) 
is  the  mean  of  Xjk  for  the  prior  data,  k  <  0 
^ik  —  is  the  ratio  of  observed  to  modeled  brightness  at  pass  k 

^mjk 

Xj  is  the  vector  of  Xj^ 

Xj^=  P(NOMj,\rj,) 

Xl  denotes  the  sample  mean  of  observations  Xji,. .  .,Xjk,,  k  >0.  If  the  Xji,. .  .,Xjk 
Xy  =  Summation  over  all  permissible  values  of  y 
z^l2  denotes  the  (100)(a/2)*  normal  percentile 

APPENDIX  B:  NUMERICAL  EXAMPLE:  CROSSTAG  DETECTION 

This  numerical  example  considers  the  scenario  for  cross-tag  detection  in  a  notional  GEO  cluster  consisting  of  a  pair 
of  cluster  peers,  Cl  and  C2.  The  Brightness  Data  for  satellites  Cl  and  C2  was  generated  using  a  numerical 
simulation  performed  by  the  Optical  Signatures  Code  [Reference  28].  The  simulation  data  emulated  one  data  point 
per  pass  where  the  mean  gap  between  the  passes  was  one  hour.  There  was  no  daytime  data.  The  simulation  data  was 
considered  for  a  period  of  two  weeks.  The  first  week  of  data  was  treated  as  prior  data  and  the  second  week  as  the 
new  data.  Two  cases  were  considered  for  the  simulations  as  shown  in  Table  Bl. 

Table  Bl:  Two  cases  of  numerical  simulation  for  CrossTag  analysis 


Case 

Satellite  Cl 

Satellite  C2 

1 

Weekl 

NOM 

NOM 

1 

Week  2 

NOM 

NOM 

2 

Weekl 

Same  as  Case  1 

Same  as  Case  1 

2 

Week  2 

CrossTag 

CrossTag 

The  prior  data  was  assumed  to  be  available  at  all  times  and  the  new  data  was  assumed  to  become  available  as  one 
Brightness  Data  point  at  a  point.  The  belief  calculation  was  performed  each  time  a  new  Brightness  Data  point  is 
received.  The  calculations  were  performed  using  the  two-facet  model  as  the  Photometry  Model,  the  Bayes  network 
and  the  control  logic  as  illustrated  in  Figure  16.1  [Reference  25].  The  Statistics  Model  was  implemented  such  that 
the  Control  Logic  can  call  any  analysis  procedure  on  an  as-needed  basis.  The  results  are  reported  in  Figures  Bl  to 
B7.  Table  B2  provides  a  key  to  how  the  figures  are  organized.  The  various  different  results  shown  in  the  figures  Bl 
to  B7  correspond  to  the  Boxes  in  Figure  16.1.  Note  that  each  underlying  calculations  is  a  discrete  step,  simple  for  a 
desktop  computer  to  execute.  Each  step  was  implemented  in  reusable  software  modules  and  the  sequence  of  steps 
was  automatically  executed  by  the  Control  Logic  [Reference  25].  The  specifics  of  each  step  performed  in  order  to 
generate  Figures  Bl  to  B7  are  given  after  Table  B2. 

\ 
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Table  B2:  Organization  key  for  the  results  in  Figures  B1  to  B7 
Week  1  is  prior  data  and  Week  2  is  new  data 


Eigure 

Satellite 

Cases 

Week 

Bla 

Cl 

1  and  2 

1 

Bib 

Cl 

1  and  2 

1 

B2a 

C2 

1  and  2 

1 

B2b 

C2 

1  and  2 

1 

B3a 

Cl 

1  and  2 

1 

B3b 

Cl 

1  and  2 

1 

B4a 

C2 

1  and  2 

1 

B4b 

C2 

1  and  2 

1 

B5a 

Cl 

1  and  2 

1 

B5b 

Cl 

1  and  2 

1 

B6a 

C2 

1  and  2 

1 

B6b 

C2 

1  and  2 

1 

B7a 

Cl 

1 

2 

B7b 

Cl 

2 

2 

B7c 

Cl 

2 

2 

Figure  Bla  shows  a  plot  of  Brightness  data  versus  time  for  satellite  Cl.  This  is  prior  data  that  spans  eight  days. 
The  last  data  point  in  this  prior  data  corresponds  to  pass  k  =  -1.  The  observed  and  modeled  Brightness  Data  is 
shown  in  Figure  Bib  for  satellite  Cl  as  function  of  the  orbit  angle.  The  modeled  data  consists  of  individual 
contributions  of  the  Body  and  the  Panels  to  the  satellite  signature.  The  orbit  angle  is  defined  in  a  manner  similar 
to  the  longitudinal  phase  angle,  except  that  the  phase  angle  is  calculated  with  respect  to  the  equatorial  plane  and 
the  orbit  angle  is  computed  with  respect  to  the  orbital  plane  of  a  satellite.  Note  that  the  Figure  Bib  is  for  the 
NOM  state  of  the  satellite.  Figure  B2a  and  B2b  show  the  corresponding  calculations  for  satellite  C2. 

The  prior  data  is  utilized  to  compute  the  prior  probability  distribution  for  Tjj^  under  NOM  conditions.  This 
calculation  is  performed  using  Equation  2.  It  corresponds  to  Box  3  in  the  Control  Logic  shown  in  Figure  16.1 
for  satellite  Cl.  The  prior  distribution  is  extracted  from  the  histogram  and  its  Gaussian  fit  shown  in  Figures  B3a 
and  B4a,  respectively  for  satellites  Cl  and  C2. 

The  ANOM  condition  is  assumed  to  only  comprise  the  CrossTag.  Or,  the  only  term  in  the  right  hand  size  of 
Equation  4  is  assumed  to  be  due  to  CrossTag.  Equation  10  is  utilized  to  determine  the  prior  probability 
distribution  under  ANOM  conditions.  This  distribution  is  also  extracted  from  the  histogram  shown  in  Eigures 
B3b  and  B4b  for  the  cluster  peers  Cl  and  C2.  The  histograms  are  quite  different  for  the  NOM  and  ANOM 
conditions  even  though  the  visual  magnitudes  for  the  two  satellites  are  quite  comparable  to  each  other.  This  is 
because  of  the  specific  character  of  the  satellite  Brightness  as  a  function  of  phase  angle  and  it  is  different  for  the 
two  peer  satellites.  The  NOM  distribution  is  similar  to  a  Gaussian  distribution,  however  it  has  longer  tails.  The 
ANOM  distribution  is  unlike  a  Gaussian  distribution.  Note  that  this  character  of  the  NOM  and  ANOM 
distribution  also  depends  on  the  accuracy  of  the  Photometry  Model,  which  varies  from  one  satellite  to  another. 
Nonetheless,  this  difference  in  the  NOM  and  ANOM  distributions  is  a  key  to  the  Statistics  Model. 

The  value  of  long  run  probability  for  NOM,  tin  is  assumed  to  be  0.90  for  each  satellite.  The  terminator- to - 
terminator  NOM-to-NOM  transition  probability  is  assumed  to  be  0.95.  The  transition  probability  matrix  for  the 
daytime  gap  is  computed  using  Equation  21.  The  pass-to-pass  transition  probability  matrix  is  computed  using 
Equation  17b.  This  is  Box  4  in  the  Control  Logic. 

Eigures  B5  considers  two  situations  for  satellite  Cl  and  Eigure  B6  considers  the  corresponding  two  situations 
for  satellite  C2.  Eigure  B5a  shows  a  plot  of  rjj^  versus  orbit  angle  for  satellite  Cl  if  it  was  correctly  tagged  in  the 
new  data  (i.e.  the  second  week  of  Brightness  Data).  Eigure  B5b  shows  a  plot  of  rjj^  versus  orbit  angle  if  the 
satellite  Cl  was  cross-tagged  with  satellite  C2  in  the  new  data.  Eigure  B6a  shows  the  plot  when  satellite  C2  was 
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correctly  tagged  in  the  new  data  and  Figure  B6b  when  it  was  cross-tagged  in  the  new  data.  The  calculations  for 
the  ANOM  histogram  are  performed  using  Equation  10. 

•  Figure  B7a  considers  the  case  when  satellite  Cl  was  correctly  tagged.  It  shows  four  graphs.  The  x-axis  is  the 
orbital  pass  number,  or  k.  The  blue  graph  is  a  plot  of  rjj^.  Note  the  peaks  in  the  values  of  rjj^  ,  which  can  be  real 
or  due  to  limitations  of  the  Photometry  Model.  The  red  line  is  the  plot  of  P(rjk  \  NOMj^  ),  which  is  obtained 
from  the  PDF  derived  from  the  histogram  shown  in  Figure  B3a.  The  green  line  is  the  plot  of  P(  \  ANOMjj^  ), 
which  is  obtained  from  the  PDF  derived  from  the  histogram  shown  in  Figure  B3b.  The  purple  line  is  the  plot  of 
P(NOMjj^  \rjk),  which  is  obtained  using  the  Bayes  theorem.  The  values  of  P(rjk  \  NOMj^  )  are  close  to  one 
for  much  of  the  time.  This  is  expected  since  the  satellite  Cl  was  correctly  tagged  in  this  case.  Towards  the  end, 
there  is  a  drop  due  to  a  consecutive  set  of  values  of  for  which  P^Tjj^  \  NOMj^  )  is  low.  This  is  Box  6  in  the 
Control  Logic  flow  chart  in  Figure  16.1.  The  persistently  high  value  of  P(  \  NOMj^  )  may  be  utilized  by  a 
user  as  indicator  for  correct  tagging  of  a  satellite. 

•  Figure  B7b  considers  the  same  calculations,  except  for  the  case  when  satellite  Cl  was  cross-tagged.  The  purple 
line  plot  for  P(  NOMj^  \  Vjj^  )  includes  large  segments  where  the  probability  of  NOM  is  close  to  zero.  This  may 
be  considered  an  indicator  by  the  user  for  the  presence  of  CrossTag  anomaly. 

•  Figure  B7c  considers  the  calculations  performed  for  the  z-score  for  P(AN0Mjj^  I  ^7c  the  cluster-based 

evidence  ,  the  model-based  evidence  and  the  belief  ^  (ANOMj^)  for  CrossTag.  This  is 

calculated  as  per  Equation  29.  The  numerical  values  and  are  not  allowed  to  become  larger  than 

a  threshold  value  (e.g.  10)  in  order  to  prevent  one  message  from  overwhelming  the  other.  The  propagated  belief 

(  ANOMj^  )  shows  a  significant  regime  for  ANOM,  which  may  be  utilized  by  a  user  as  indicator  for  the  of 
ANOM  conditions.  This  is  Box  8  in  Figure  16.1.  Figure  B7c  shows  a  notional  red  line,  which  is  meant  to  denote 
a  user  defined  threshold  for  the  presence  of  CrossTag. 

•  Note  that  the  (  ANOMj^  )  drop  significantly  after  several  passes  (Figure  B7c).  This  is  in  conflict  with  the  high 
values  of  belief  computed  for  the  balance  of  the  passes.  This  occurs  because  the  Cross -tagged  satellites  are 
equally  bright,  which  is  makes  it  difficult  to  distinguish  between  them  on  the  basis  of  their  Brightness. 

•  If  the  level  of  belief  for  NOM  in  the  new  data  is  on  the  same  level  as  the  prior  data  and  it  exceeds  a  user -defined 
threshold,  the  time  slider  is  shifted  forward  to  pass  k.  This  is  Box  9  in  the  Control  Logic. 

Note  that  Figures  B1  to  B6  represent  the  calculations  performed  on  the  prior  data.  These  calculations  need  to  be 
performed  each  time  the  time  slider  is  advanced  (Section  15).  It  takes  several  minutes  of  computations  in  order  to 
perform  these  calculations.  Majority  of  this  time  is  consumed  by  the  Inversion  Model.  Figure  B7  calculations  would 
be  performed  each  time  a  new  point  of  Brightness  Data  is  received.  These  calculations  take  only  a  few  seconds  to 
complete  and  thus  the  evolving  belief  in  NOM  or  ANOM  status  can  be  obtained  in  near  real  time. 

This  example  considered  an  idealized  scenario  where  the  cluster  comprised  only  two  satellites  and  the  only  anomaly 
could  occur  was  CrossTag.  In  a  more  realistic  scenario,  the  cluster  may  contain  more  satellites  and  all  anomalies  are 
possible.  The  calculations  necessary  to  analyze  such  a  scenario  consist  of  many  more  of  that  shown  in  Figure  B7. 
Figures  B1  to  B6  are  for  the  prior  data  and  they  are  common  to  all  ANOM  assessments. 

A  salient  next  step  is  to  synthesize  the  evolving  belief  automatically  in  order  to  rank  the  potential  ANOM  situation. 
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Figure  B la:  Brightness  Data  vs.  Time  for  satellite  Cl  (prior  data) 
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Figure  Bib:  Observed  and  Predicted  (Modeled)  Brightness  Data  for  Satellite  Cl 
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Figure  B2a:  Brightness  Data  vs.  Time  for  satellite  C2  (prior  data) 


Figure  B2b:  Observed  and  Predicted  (Modeled)  Brightness  Data  for  Satellite  C2 
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Figure  B3a:  Prior  distribution  /(rjk  |  NOMf)  for  satellite  Cl 
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Figure  B3b:  Prior  distribution  f{rju  \  ANOMj)  for  satellite  Cl 
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Figure  B4a:  Prior  distribution  f{rjk  \  NOMj)  for  satellite  C2 
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Figure  B4b:  Prior  distribution  f(rjj^  \  ANOMj)  for  satellite  C2 
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Figure  B5a:  Plot  of  Vjj^  for  new  data  (k  >  0)  versus  orbit  angle  for  satellite  Cl  (NOM) 
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Figure  B5b:  Plot  of  for  new  data  (k  >  0)  versus  orbit  angle  for  satellite  Cl  (ANOM) 
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Figure  B6a:  Plot  of  for  new  data  (k  >  0)  versus  orbit  angle  for  satellite  C2  (NOM) 
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Figure  B6b:  Plot  of  rjj^  for  new  data  (k  >  0)  versus  orbit  angle  for  satellite  C2  (ANOM) 
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Figure  B7a:  Plot  of  P{NOMj^  \  )  for  new  data  (k  >  0)  when  satellite  Cl  is  correctly  tagged 
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Figure  B7b:  Plot  of  P(N0Mjj^  \  Tjj^  )  for  new  data  (k  >  0)  when  satellite  Cl  is  cross-tagged 
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